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ABSTRACT 


The  mathematical  pendulum  is  one  of  the  most  widely  studied  problems  in 
engineering  physics.  However,  this  is  primarily  limited  to  the  classical  pendulum 
with  a  single  massless  bar  and  mass  configuration.  Extensions  to  this  include 
multi-degree  of  freedom  systems,  but  many  of  the  classical  assumptions,  such  as 
a  single  bar  per  mass,  are  preserved.  Several  designs  used  in  practice  utilize 
multiple  bars  or  trapezoidal  configurations  to  enhance  stability.  Such  designs 
have  not  been  studied  in  detail.  Also,  there  is  a  need  for  additional  work  to  fully 
analyze  their  response  characteristics.  The  two-string  pendulum  design 
characteristics  are  initially  investigated,  both  in  terms  of  oscillation  and  string 
tension.  Analytical  and  numerical  methodologies  are  applied  to  predict  the 
response  of  the  two-string  pendulum  in  free  and  forced  oscillations.  Results  are 
validated  utilizing  comparisons  generated  by  simulations  conducted  with  a 
standard  commercial  software  package.  A  preliminary  optimization  study  is 
conducted  for  a  driven  two-string  pendulum.  Finally,  in  a  typical  design  case,  it  is 
shown  how  to  apply  the  results  of  the  analysis  and  optimization  studies 
developed  in  this  work. 
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I. 


INTRODUCTION 


A.  BACKGROUND 

The  pendulum  is  one  of  the  most  commonly  studied  mechanical 
structures.  Galileo  himself  made  the  first  observation  of  the  pendulum’s 
performance.  He  set  the  roots  of  the  pendulum’s  physics  by  realizing  that  the 
frequency  of  oscillation  for  a  given  pendulum  is  independent  of  its  displacement’s 
amplitude  for  small  amplitude  oscillations,  which  enables  the  linearization  of  the 
equations  of  motion.  [1], 

Following  this  initial  observation,  numerous  studies  were  performed  on  the 
pendulum’s  physics.  These  studies  were  not  limited  to  the  simple  or  standard 
pendulum,  but  were  extended  to  more  complex  pendulum  structures.  Some 
examples  are  the  bifilar  pendulum,  the  torsion  pendulum,  the  double  pendulum, 
and  the  Foucault  pendulum. 

To  describe  analytically  the  oscillatory  behavior  of  the  numerous  types  of 
pendula,  several  techniques  were  developed.  These  techniques  attempted  to 
predict  the  solution  of  the  pendulum’s  dynamics,  either  by  treating  its  attitude  as 
a  linear  approximation,  which  produced  acceptable  and  fairly  accurate  results  for 
a  limited  number  of  cases  (linearized  pendulum,  ,  etc.),  or  by  facing  its  nonlinear 
character  recruiting  and  applying  complex  methods  (perturbation  methods  or  the 
Lagrange  method,  which  produces  the  exact  nonlinear  equation  of  motion). 

From  a  mathematical  analysis  viewpoint,  the  more  composite  the  structure 
of  the  pendulum,  the  more  complex  the  solution  approach  method.  Some  types 
of  pendula  present  a  linear  character  when  the  oscillation  is  constrained  in  small 
displacement  amplitudes  (standard  and  double  pendulum).  Large  amplitude 
oscillations  are  nonlinear  and  sometimes  show  chaotic  behavior  (double 
pendulum).  If  the  problem  is  approached  from  an  energy  perspective,  in  linear  or 
linearized  undriven  oscillations,  energy  is  conserved.  Nevertheless,  there  are 
examples  where  energy  is  not  conserved.  Such  cases  may  be  simple  to  analyze, 
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such  as  an  oscillation  with  a  constant  damping  coefficient.  Others,  such  as 
oscillations,  including  jumps  or  variable  accelerating  support  motion,  are  not 
simple  to  analyze. 

B.  MOTIVATION  -  OBJECTIVE 

Several  mechanical  structures  have  dynamic  behavior  similar  to  that  of  a 
pendulum.  Some  examples  are  different  types  of  cranes  and  crane-like 
mechanisms  used  for  marine  vehicles,  such  as  AUVs,  rescue  boats,  and  naval 
vessels  for  launching  and  recovery.  Their  multi-rope  suspension  mechanism  is 
not  designed  to  oscillate;  rather,  it  is  designed  to  prevent  swinging  and  to 
promote  stability.  Another  example  of  a  pendulum  system  designed  to  prevent 
swinging  is  the  Blackburn  pendulum.  [1], 


A 


Figure  1.  The  Blackburn  pendulum 

The  Blackburn  pendulum’s  suspension  system  is  comprised  of  two  strings 
forming  a  "Y”  shape.  The  purpose  of  this  design  is  to  provide  two  distinct  and 
effective  pendulum  string  lengths.  Therefore,  for  an  oscillatory  motion  within  a 
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plane  defined  by  the  Y  shape,  the  length  of  the  pendulum  is  the  distance  from 
point  B  to  point  C.  In  contrast,  for  a  motion  within  a  plane  perpendicular  to  the  Y 
shape,  the  effective  length  is  the  distance  from  point  A  to  point  C.  Therefore,  two 
distinct  frequencies  of  oscillation  define  the  motion  of  the  Blackburn  pendulum  in 
the  two  perpendicular  planes.  Thus,  its  modes  of  oscillation  are  described  by 
superposition  of  the  two  basic  frequencies. 

Nevertheless,  all  the  previously  mentioned  characteristics  of  the  Blackburn 
pendulum  depend  on  the  assumption  that  point  B  will  not  be  able  to  oscillate  with 
respect  to  one  of  the  two  strings  that  starts  from  B  and  ends  at  the  support 
points,  i.e.,  the  distance  of  B  from  the  two  support  points  will  always  be  constant. 
However,  this  assumption  is  valid  under  restrictions.  A  violent  motion  of  the 
support  would  theoretically  be  able  to  disturb  the  partial  equilibrium  of  the  "Y” 
pendulum. 

The  motivation  for  this  thesis  was  to  obtain  analytical  expressions,  as  well 
as  numerical  solutions,  to  describe  the  behavior  of  oscillating  pendula  with 
designs  similar  to  the  one  of  the  "Y”  pendulum.  The  basic  model  design  that  will 
be  initially  analyzed  is  the  one  of  a  Blackburn  pendulum  in  which  point  A 
coincides  with  B,  and  is  allowed  to  oscillate  only  within  the  plane  defined  by  this 
“V”  shape  (more  detailed  description  will  be  provided  in  the  following  sections). 
Moreover,  since  such  a  "V”  design  is  the  basic  element  component  of  more 
complex  suspension  designs,  its  capability  of  preventing  or  damping  oscillations 
will  be  investigated. 
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THE  TWO-STRING  PENDULUM 


A.  MATHEMATICAL  DESCRIPTION  OF  THE  TWO-STRING  PENDULUM 

1.  Introduction 

The  purpose  of  this  section  is  to  analyze  the  motion  equations  of  the  two- 
string  pendulum.  For  brevity  within  this  thesis,  this  will  be  referred  to  as  the  2- 
pendulum.  Several  techniques  will  be  used  to  analyze  the  motion  of  the  2- 
pendulum  with  the  aim  to  present  an  accurate  solution  to  the  nonlinear  equations 
of  motion  that  describes  its  motion.  The  main  purpose  is  to  compare  main 
features  of  the  dynamic  response  of  the  2-pendulum  with  those  of  the  standard 
single-mass,  single-string  pendulum.  These  comparisons  will  be  in  terms  of  the 
amplitude  of  the  resulting  motion,  the  resulting  natural  frequency  of  the  oscillator, 
as  well  as  the  string  tension.  Different  values  of  initial  velocity  and  initial  angular 
displacement  will  be  used  to  investigate  the  system’s  response  to  both  an  initial 
excitation  and  an  initial  displacement  from  its  equilibrium  position. 


2.  The  Standard  Pendulum 


Figure  2.  Small  angle  oscillation  of  the  standard  pendulum 
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First,  we  will  investigate  the  motion  of  the  simple  gravity  pendulum 
undergoing  free  undamped  vibration.  The  model  that  we  will  use  is  shown  in 
Figure  1  —  a  sphere  of  mass  m  which  is  suspended  with  a  massless  string  of 
length  =  i  from  its  center  of  gravity. 

To  derive  the  equations  of  motion,  we  will  use  the  moment  method. 
According  to  Newton’s  Law,  the  inertial  moment  of  the  mass  m  must  be  equal  to 
the  sum  of  moments  applied  on  the  mass. 

mi2  6  =  -mg£  sin  0  (2.1) 

mPG  +  mgl  sin  0  =  0  (2.2) 

Equation  (2)  above  is  nonlinear  since  it  contains  the  term  sin#.  We  can 
substitute  the  trigonometric  nonlinearity  with  its  Taylor  series  expansion 

n3  n5 

sin0  =  0-|^+|p..  (2.3) 

To  simplify  our  calculations,  we  will  initially  assume  only  small  amplitude 
oscillations.  In  this  case,  we  can  assume  without  compromising  accuracy  that 

sin  #  «  #  (2.4) 
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Figure  3.  Linear-nonlinear  discrepancy 

As  demonstrated  by  Figure  3,  there  is  a  range  of  Q  where  equation  (4)  is 
valid  and  provides  a  good  approximation.  As  depicted  in  the  same  figure,  the 
difference  between  the  linearized  and  the  nonlinear  expression  is  negligible  for 

TC 

values  of  9  <  — . 

10 

For  such  small  amplitude  oscillations,  the  equation  of  motion  can  be 
linearized 


m£29  +  mg£9  -  0 

(2.5) 

o 

ll 

+ 

(2.6) 

The  solution  to  the  above  differential  equation  for  initial  displacement  and 
velocity  of 
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becomes 


0(0)  =  0O 
0(0) =4 


(2.7) 


(9  = 


M 


sin  cont  +  0(0)  cos  cont 


(2.8) 


where 


(2.9) 


is  the  natural  frequency  of  the  linear  pendulum. 
From  Figure  2,  the  restoring  moment  is 

Fr  =  mg  sin  0 


or,  for  small  oscillations, 


Fr  =  mg6 


(2.10) 


(2.11) 


Another  feature  of  the  pendulum  worth  mentioning  is  the  pendulum  jump 
[5],  Pendulum  jumps  occur  when  the  mass  leaves  the  circular  path  defined  by  the 
string  and  follows  other  trajectories.  Such  response  is  only  observed  if  the  string 
tension  becomes  zero.  This  is  dependent  on  the  position  of  the  pendulum  and 
the  initial  conditions.  From  Figure  2,  we  can  visualize  that  the  centripetal  force, 
which  obliges  the  mass  to  follow  a  circular  path,  is 

2  2 

T- t  THU  rri  /)  „  „  rji  THU  /j  /,-j  A  \ 

Fc  =  — - —  ~T  —  mg  cos  6  <^>  T  =  — - — + mg  cos  $  (2.12) 


where  u  is  the  linear  velocity  of  the  mass.  The  magnitude  of  u  depends  only  on 
the  initial  conditions,  since  only  conservative  forces  govern  the  pendulum’s 
motion.  Assuming  the  initial  conditions  (2.7),  which  are  equivalent  to 

zQ  =  £(\  —  cos  0q)  and  uo=10o,  the  initial  energy  of  the  mass  is 


Eo=To+Vo=\mul 


mgzQ.  During  oscillation,  the  energy  of  the  mass  for  an 
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|  r\ 

arbitrary  angle  G  would  be  E  =  T  +  V  =  —mu  +  mgz  ,  where  z  =  £(\-cos0) . 

Since  the  total  energy  is  conserved  throughout  the  oscillation  (absence  of 
damping  or  other  non-conservative  forces)  E  =  E0 ,  the  magnitude  of  the  linear 

velocity  would  be 

E 

n  =2.  gz  (2.13) 

V  m 


Finally,  the  tension  of  the  string  throughout  the  oscillation,  from  (2.12)  and 
(2.13),  will  be  equal  to 


2(E0  -  mgz) 

T  = - — ^—+mg  cos  6 


(2.14) 


where  a  positive  sign  denotes  orientation  towards  the  center  of  oscillation.  This 
promotes  the  circular  path  motion  of  the  mass.  The  only  term  of  the  latest 
equation  that  may  take  negative  values  and,  thus,  lead  to  zeroing  the  string 
tension,  is  mg  cos  0 .  It  may  only  happen  for  angular  displacements  greater  than 


n 

2  ' 

The  possibility  of  jumps  is  present  only  if  the  initial  conditions  allow  the 

oscillator  to  climb  to  0>  — ■  After  that  point,  the  path  of  the  mass  would  be 

governed  by  the  regular  equations  which  define  a  projectile  with  a  given  initial 
linear  velocity  magnitude  and  orientation;  that  is,  until  the  position  of  the  mass 
will  again  reach  a  distance  equal  to  the  natural  length  of  the  pendulum’s  string. 
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3. 


The  2-Pendulum 


Figure  4.  The  2-pendulum 


The  geometry  of  the  2-pendulum  is  pictured  in  Figure  4.  The  main 
distinguishable  characteristic  of  the  2-pendulum,  compared  to  the  simple 
pendulum,  is  that  the  mass  m  is  now  suspended  by  two  strings.  To  be  consistent 
with  the  previous  model,  the  strings’  perpendicular  distance  from  the  plane  of 

£ 

suspension  will  be  / .  Consequently,  /  = - .  The  point  of  application  of  the 

cos  a 

strings  on  the  mass  is  again  at  its  center  of  gravity  while  each  string  forms  an 
angle  a  with  respect  to  the  perpendicular.  Our  aim  is  to  evaluate  the  effect  of  this 
angle  on  the  motion  of  the  mass  m. 

For  the  purposes  of  this  preliminary  analysis,  we  will  assume  that  the 
governing  force  during  the  pendulum  oscillation  is  the  weight.  That  means  that 
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we  will  not,  for  now  examine  the  forces  applied  on  the  mass  to  the  point  where  it 
reaches  the  lower  point  of  its  trajectory  and  switches  the  string  about  which  it 
swings. 


gcos(0+a 


Figure  5.  Small  angle  oscillation  of  the  2-pendulum,  0<O 
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Assuming  an  angle  9  (counterclockwise  for  positive  values),  as  in  Figure 
6,  and  equating  inertia  and  non-inertia  forces  applied  on  the  mass  as  before,  the 
resulting  equation  of  motion  is 

ml2ij  +  mglsm(0  +  a^  =  0  (2.15) 

or 

#  +  ysin(#  +  a)  =  0  (2.16) 

To  further  develop  equation  (2.16),  we  should  take  into  account  that  the 
term  sin(#  +  a)  takes  positive  values  for  the  case  of  0>O  (Figure  6),  and 

negative  values  for  the  case  of  9< 0  (Figure  5).  However,  since  the  restoring  force 
should  always  oppose  further  increase  in  the  magnitude  of  the  angle  9,  it  takes 
the  same  magnitude  for  the  same||a  +  0||.  Furthermore,  using  the  trigonometric 

identity  sin(a  +  6)  =  sin <2 cos 6  + cos <2 sin 6 ,  equation  (2.16)  becomes 


#  +  y(sin6>coscr  +  cos#sincr)  =  0 


(2.17) 
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Examining  equation  (2.17),  we  can  see  that  the  consistency  of  the 
restoring  force  orientation  is  preserved  in  the  term  sin  looser  (because  of 
sin#),  but  not  in  sin  a  cos#  (since  cos#  does  not  change  sign  when  G  changes 
sign).  Consequently,  to  preserve  the  restoring  force  consistency  for  all  values  of 
Q,  we  should  rewrite  equation  (2.17)  as  follows 


#  +  ySin#cos<2±yCOS#sina  =  0 


This  is  equivalent  to 


#  +  ysin#cosa  +  yCOS#sin«sign#  =  0 


g. 


(2.18) 


The  sign  nonlinearity  assumes  the  value  (+1)  for  positive  and  (-1)  for 
negative  arguments.  In  addition  to  the  above,  if  we  assume  small  angles  of  6,  we 
can  rewrite  (2.18)  as  follows 

#  +  y#costt  +  ySin<2sign#  =  0  (2.19) 

The  restoring  force  from  equation  (2.15)  is 

Fr  =  mg  sin(#  +  a)  (2.20) 

For  small  oscillations,  equation  (2.20)  becomes 

Fr  =  mgO  cos  a  +  mg  sin  a  sign#  (2.21) 

Our  goal  is  to  solve  analytically  and  numerically  differential  equations 

(2.18). 

Finally,  based  on  the  previous  analysis  performed  on  the  jumps  for  the 
standard  pendulum,  we  realize  that  the  2-pendulum  may  exhibit  the  same 

71  71 

behavior,  only  now  the  required  angular  displacement  is#  +  er  >y  <=>  9>  —  -a  . 


4.  Comments 

Despite  its  similarity  with  the  standard  pendulum,  equation  (2.18)  exhibits 
some  major  characteristic  differences.  The  most  significant  difference  between 
the  standard  pendulum  and  the  2-pendulum  is  that  the  latter  possesses  a 


13 


discontinuity  in  its  stiffness  around  the  equilibrium  position.  As  can  be  seen  from 
(2.19),  there  is  an  instantaneous  jump  in  the  value  of  the  restoring  moment 
around  zero.  This  is  depicted  in  Figure  7  where  the  normalized  spring  constant  of 
the  2-pendulum  is  shown  for  different  values  of  a ,  along  with  the  linear  spring 
constant  of  the  standard  pendulum.  This  discontinuity  at  zero  prevents  us  from  a 
classical  linearization  of  the  2-pendulum  and  requires  different  techniques  for  the 
analysis.  Such  techniques,  based  on  the  perturbation  theory,  are  the  subject  of 
the  analysis  that  follows. 

1.5 

1 

0.5 

0 

-0.5 

-1 

-1.5 

-1  -0.8  -0.6  -0.4  -0.2  0  0.2  0.4  0.6  0.8  1 

Figure  7.  Normalized  spring  constant  of  the  2-pendulum 


Normalized  Spring  Constant 


i - i - 1 - \ - 1 - I - ! - 1 - r 


J _ I _ I _ I _ I _ I _ I _ I _ L 


14 


5. 


Velocity  and  Force  Vector  Analysis 


Figure  8.  General  force  and  velocity  diagram  at  the  switch  point 

Before  advancing  with  the  solution  of  equations  (2.18)  and  (2.19),  we 
should  take  into  account  the  dynamics  of  the  mass  m  at  the  point  when  it 
switches  strings  about  which  it  swings.  In  Figure  8,  we  can  visualize  the  forces, 
weight  B,  instantaneous  string  tension  T,nst  (blue  with  continuous  line  vectors), 
and  their  resultant  R  (blue  dashed  line  vector).  During  the  change,  it  acts  on  the 
mass  and  actually  contributes  to  not  only  the  switch  at  the  center  of  rotation,  but 

to  the  velocities  just  before  Vb  (red  vector),  and  just  after,  Va  (green  vector)  -- 
the  switch.  Nevertheless,  since  the  weight  of  the  mass  is  a  conservative  force 
and  does  not  undergo  changes  in  magnitude  during  this  switch,  we  can  safely 
assume  that  its  contribution  in  the  velocity  direction  switch  is  negligible. 


15 


Figure  9. 


Force  and  velocity  diagram  at  the  switch  point  for  o= 45° 


According  to  the  vector  diagram  in  Figure  9,  there  is  a  limiting  value  of 
angle  a,  which,  if  exceeded,  the  string  switch  is  impossible.  In  the  instance  that 

Tjnst  lies  in  the  same  direction,  but  in  opposite  orientation  to  the  Vb ,  this  force  is 


unable  to  change  the  direction  of  the  velocity  from  vb  to  Va ,  irrespective  of  its 
magnitude.  This  limiting  value  is  shown  to  be  a=45°.  In  this  case,  this  exact  value 
of  the  mass  will  either  stop  or  bounce  in  the  same  direction.  However,  the 


orientation  of  Tjnst,  will  have  a  velocity  magnitude  less  than  or  equal  to  Vb 
(depending  on  the  energy  absorption  of  the  string).  As  a  result,  the  following 
preliminary  analysis,  Chapter  II,  Section  A,  numbers  6,  7,  and  8,  will  be  valid  only 
for  cr<45°.  Furthermore,  we  will  not  consider  the  exact  forces  that  affect  dynamics 
and  energy  equivalence  at  the  switch  point.  Rather,  we  will  only  take  into  account 


that  their  effect  is  the  jump  of  the  velocity  from  Vb  to  Va,  where 


Va 


6.  Numerical  Integration  Solution 

The  numerical  solution  is  performed  using  the  Matlab  function  ode4  5. 
This  ode4  5  Matlab  function  solves  numerically  the  equation  of  motion  and  is  able 
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to  take  into  account  nonlinearities,  such  as  trigonometric  functions  and  sign 
changes.  Thus,  ode4  5  will  provide  the  exact  solution,  provided  that  an 
appropriate,  fine  enough,  time  step  is  selected.  The  PrograrrM.m  Matlab-file 
was  used  to  plot  the  angular  displacement  6  versus  time  for  different  values  of 

angle  a,  assuming  initial  conditions  0(0)  =  0O  and  0(0)  =  0.  In  parallel,  we 

plotted  the  angular  displacement  6  versus  time  for  a  simple  pendulum  of  string 
length  i .  In  addition,  the  Program_2.m  provided  the  same  kind  of  plots,  where 
the  initial  conditions  were  0(0)  =  0  and  0(0)  =  0O  . 


7.  Analytical  Solution 

Three  different  methods  will  be  used  to  attempt  to  derive  the  exact  or 
approximate  solution  to  equation  (14)  or  its  equivalent  (16). 


a.  Exact  Solution  for  a  Time  Increment 

One  method  we  could  use  to  obtain  the  solution  of  equation  (2.18) 

is  to  solve  it  for  the  time  period  where  the  term  y  sin  a  sign  0  has  a  constant 

sign  [2],  Thus,  assuming  an  initial  angular  displacement  @(0)=A,  the  equation  of 
motion  for  0<@<A  is 


or 


g  . 

sm  a 

0  =  --t - + 


r 


g_ 

l 


cos  a 


A  + 


g .  ^ 

ysm  a 

g 

^-cos  a 
/  7 


cos 


|yC0S<2 


"A 


(2.22) 


0  =  -  tan  a  +  (A  +  tan  a)  cos 


Mucosa 


/ 


t 


(2.23) 


Let  us  name  the  square  of  the  pseudo-natural  frequency  of 
oscillation  and  of  the  equivalent  pseudo-forcing  amplitude  as  bellow 


2  8 

COq  =yC0S<2 

P-^-sina 


(2.24) 
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At  the  time  t  =  —  arccos 


co. 


o 


f  \ 

1 

1 1  M 

p  j 


when  0=0,  the  velocity  is 


u(tx)  =  -coQ 

and  then  the  equation  of  motion  becomes 


r  ,  P  A 
A  +  — y 

V 


sm  c o0t 


with  initial  conditions  ii{tl)  =  -coQ 


6  +  q)qO  =  P 

r  A 

A  H - y 

\  )J 


(2.25) 

sin co0t  and  0(tx)  =  0,  and  so  on.  The 


value  t{  is  equivalent  to  one  quarter  of  the  period  of  oscillation. 
Finally  the  period  of  oscillation  is 


T  = — arccos 
exact  oo 


f  \ 

1 

1 1  M 

p  j 


or 


exact 


g 

ycos  a 
t 


-  arccos 


f  \ 

1 

1  +  —^— 
v  tan  a  ) 


so  the  frequency  of  oscillation  is 


co 


2  71 


exact  T 


exact 


(2.26) 


(2.27) 


(2.28) 


b.  Method  of  Slowly  Changing  Phase  and  Amplitude 

A  second  method  for  obtaining  the  analytic  solution  only  for  the 
linearized  equation  (2.19)  or ,  considering  (2.24) 

O  +  colO+P  sign6>  =  0,  (2.29) 
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is  using  the  method  of  Slowly  Changing  Phase  and  Amplitude  [2],  This  method 
states  that  the  equations  of  the  slowly  changing  phase  y/  and  amplitude  a  are 


=  *ffl,(,)Jo2j[-/’{asin(6'  +  ^)}]sin(a  +  y)^ 

i  =  — — j — r-\^^-Psign[asin[0+i//)^sm(a+i//)d2 


(2.30) 

(2.31) 


Evaluating  the  above  integrals, 


a  =  0  and  = 


4  P 


IndiCOr, 


a  =  A  and  y/  = 


IP 

7lA(Dr 


■  +  ^o 


;r 


Finally,  if  we  choose  the  initial  conditions,  such  as  —  ,  the  approximate 


solution  is 

6  =  Acos 

and  the  approximate  period  is 

T . =  - 


2  P 

(>\)  + - 7 — 

yr/ky 


oy 


2;r 


app 


ox 


A 


1  + 


2P 

TTor  A 

o  y 


(2.32) 


(2.33) 


while  the  approximate  angular  frequency  is 


® app  % 


1  + 


2  P 

7TO)qA  j 


considering  (2.24) 


a>app=\  ycosa 


r 


2  tana 


v1  +  _  xA 


(2.34) 


(2.35) 


c.  Energy  Method 

The  application  of  this  method  suggests  that  we  could  approximate 
the  equation  of  motion  of  the  enhanced  pendulum  with  one  of  a  simple 

pendulum.  However,  the  restoring  force  would  be  increased  by  a  factor  so  that 
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the  total  work  produced  by  this  force  along  a  time  period  would  be  the  same  as 
the  forces  that  govern  the  motion  of  the  enhanced  pendulum. 

From  equation  (2.20),  the  work  produced  by  the  restoring  force  of 
the  enhanced  pendulum  fora  motion  0 <O<60  is  equal  to 

'0mgsin(a  +  6)d6  (2.36) 

while  the  one  of  the  simple  pendulum  is 

mg  sin  6d 6  (2.37) 

Dividing  (2.36)  by  (2.37),  we  obtain  a  factor  of 

F  =  \°0mgsin{a  +  6)dO  /  j^0  mgsvnOdO  (2.38) 

which  provides  us  with  an  order  of  increased  magnitude  of  forces  applied  during 
the  oscillation  of  the  enhanced  pendulum  compared  to  the  ones  applied  on  the 
simple  pendulum. 

Thus,  we  will  multiply  the  term  mg  sin  6  by  F  and  solve  the 
approximate  equation  of  motion 

0  +  F^-  sin6>  =  0  (2.39) 

or  in  linear  form, 

0  +  F^G  =  0  (2.40) 

Thus,  the  approximate  angular  frequency  of  oscillation  is 


co. 


app 


(2.41) 


8.  Evaluation  of  Results 

To  evaluate  both  numerical  and  analytical  approximations,  we  will  first  plot 
the  three  angular  displacement  solutions  (ode45  numerical  integration  -  Slowly 
Changing  Phase  and  Amplitude  method  -  Energy  method).  This  will  be 
compared  to  the  solution  of  the  standard  pendulum  in  terms  of  gain  and 
frequency  of  oscillation.  This  will  best  show  that  the  2-pendulum  offers  significant 
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gains  in  terms  of  increased  restoring  moment  and  frequency  of  oscillation.  It, 
also,  measures  how  close  it  comes  to  the  exact  solution  (represented  by  the 
ode45  solution).  The  predicted  response  will  be  presented  for  two  distinct  initial 
conditions  ~  initial  angular  displacement  and  initial  angular  velocity.  In  addition  to 
the  above,  we  will  plot  the  phase  trajectories  for  the  same  initial  conditions,  as 
well  as  the  comparison  of  the  predicted  angular  frequencies  of  oscillation.  This 
will  help  to  better  visualize  the  discrepancy  of  the  prediction  methods. 
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Figure  10.  a=Tr/6 
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Figure  11.  a=Tr/12 


Figure  12.  a=Tr/18 
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Figure  13.  a=Tr/36 
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Figure  14.  a=Tr/90 
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Figure  16.  Phase  diagram  a=Tr/6 
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Figure  17.  Phase  diagram  ar=Tr/12 
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Figure  18.  Phase  diagram  a=nY18 
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Figure  19.  Phase  diagram  a=Tr/36 
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Figure  20.  Phase  diagram  a=Tr/90 
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Figure  21 .  Phase  diagram  cr=Tr/180 


Figure  22.  o= tt/6  with  initial  velocity 
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Figure  23.  a=Tr/12  with  initial  velocity 
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Figure  24.  a=Tr/18  with  initial  velocity 


28 


9  (rad)  0  (rad) 


ODE45 

—  SCPA  method 

— 

Energy  method 

Standard  pendulum 

j_ i_ i_ i_ i_ i 


0  0.5  1  1.5  2  2.5  3 

t  (sec) 

Figure  26.  a=Tr/90  with  initial  velocity 
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Figure  27.  cr=Tr/180  with  initial  velocity 
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Figure  28. 
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Figure  29.  0o=tt/1O 
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Figure  30.  @0=Tr/20 
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In  Figures  28,  29,  30,  and  31,  we  can  evaluate  the  accuracy  of  the 
response  prediction  method  with  respect  to  the  initial  displacement  6o  and  value 
of  the  characteristic  angle  o.  According  to  equations  (2.28),  (2.35),  and  (2.41), 
the  oscillation  angular  frequency  depends  on  both  parameters.  Therefore,  the 
four  plots  consider  four  arbitrary  values  for  the  initial  angular  displacements.  This 
allows  us  to  better  appreciate  the  accuracy  of  the  two  prediction  methods,  SCPA 
and  Energy,  with  respect  to  the  main  attribute  of  the  2-pendulum,  which  is  angle 
a.  We  should  note  that  the  three  equations  are  valid  only  for  the  linear 
approximations  of  the  equations  of  motion.  In  addition,  the  graphs  consider 
cr<45°,  a  restriction  we  already  mentioned  in  the  beginning  of  the  analysis. 

According  to  the  plots,  we  conclude  the  following: 
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•  The  angular  frequency  go  of  the  2-pendulum  is  greater  compared  to 
the  frequency  of  the  standard  pendulum  which  has  the  same  equivalent  string 
length  £  . 

•  A  smaller  maximum  angular  displacement  for  the  2-pendulum 
compared  to  the  standard  pendulum  results  from  applying  the  same  initial 
angular  velocity. 

•  The  discrepancy  between  the  two  predictions,  and  the  exact 
calculations,  is  very  low  for  small  values  of  a,  but  increases  for  larger  ones.  More 
specifically,  while  the  Energy  method  follows  the  exact  solution  slope  trend 
(although  providing  higher  predicted  values),  the  SCPA  method  increases  almost 
linearly. 

•  The  SCPA  method  tends  to  provide  more  accurate  predictions  for 
larger  values  of  Go.  If  Go  becomes  too  small,  the  discrepancy  from  the  exact 
solution  grows  and,  therefore,  the  predictions  are  incorrect. 

9.  Horizontal  and  Vertical  Displacements 

Despite  the  obvious  advantage  that  the  2-pendulum  possesses  in  terms  of 
increased  angular  frequency  of  oscillation,  increased  restoring  force,  and 
reduced  maximum  angular  displacements,  it  is  worth  mentioning  a  disadvantage: 
the  differences  of  both  horizontal  and  vertical  displacements  associated  with 
angular  displacements.  The  following  figure  depicts  the  relationship  that  governs 
the  maximum  horizontal  and  vertical  displacements  for  Qo= tt/18  and  a  range  of  a. 
This  figure  depicts  that,  for  each  angular  displacement  of  the  2-pendulum,  the 
horizontal  and  vertical  displacements  are  smaller  and  greater  respectively  when 
compared  to  the  standard  pendulum.  This  is  due  to  the  two  pendula’s  different 
reference  systems  used  to  measure  their  angular  displacements.  Therefore,  to 
obtain  a  better  appreciation  of  the  comparison  between  the  standard  and  the  2- 
pendulum,  the  previous  examples  of  a  will  be  plotted  in  terms  of  linear 
displacements  instead  of  angular  displacements  (Figure  33). 
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Figure  32.  Variation  of  maximum  vertical-horizontal  displacements  wrt 

characteristic  angle  a 
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Figure  33.  2-pendulum’s  mass  trajectories  during  oscillation 
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B.  STRING  TENSION  INVESTIGATION 


1.  Analytical  Formulation 

Another  important  concept  that  requires  further  investigation  is  the  tension 
applied  on  the  string  of  the  2-pendulum.  As  already  seen  in  equations  (2.15)  and 
(2.2),  the  restoring  moment  on  the  mass  m  is  larger  compared  to  the  standard 
pendulum.  Contrarily,  the  string  tension  has  smaller  values  since  two  strings, 
instead  of  one,  are  used  to  suspend  the  mass.  As  pictured  in  Figure  3,  while  in 
equilibrium,  positioning  the  string  tension  is  dependant  only  on  the  value  of  the 
characteristic  angle  a 


Ts  = 


mg 

2  cos  a 


(2.42) 


According  to  Figure  4,  during  swinging  and  while  the  sign  of  angle  6  is 
continuous,  the  string  tension  can  be  computed  as  follows 


Ts  =  Fc  +  \\mg  cos  (a +  6)j  (2.43) 

2 

TYIVI  *  9 

where  Fc  =  — - —  =  m02l  is  the  centripetal  force. 

Equation  (2.43)  shows  that  the  magnitude  of  static  Ts  ranges  from 

TVlg  n  71 

< 7s <oo,  for  0 <a  +  0 <  —  ,  proving  that  increasing  the  characteristic  angle 

does  not  have  only  positive  results  in  the  overall  design. 

To  determine  the  dynamic  string  tension,  trigonometric  equalities  are 
applied  and 

Ts  =  m02l  + 1| mg  cos  a  cos  6  -  mg  sin  a  sin  6\ 

7T 

To  attempt  to  linearize  the  above  equation  for  Q<  — , 


Ts  =  md2l  +  mg  cos  asignO  -  mg  sin  a 6 


(2.44) 
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For  the  standard  pendulum,  the  calculations  are  easier.  Therefore,  the 
tension  applied  on  the  string,  while  in  equilibrium  position,  is 


Ts  =  mg 

(2.45) 

and  during  the  oscillation  is 

Ts  =  m02l  +  \mg  cos  6\ 

(2.46) 

or,  for  small  amplitude  oscillations, 

Ts  =  m02l  +  mg 

(2.47) 

Since  ode45  solves  the  exact  differential  equation  of  motion  and,  we 
already  evaluated  the  accuracy  of  the  predictions  provided  by  the  SCPA  and 
Energy  methods,  the  string  tension  simulation  will  be  performed  using  the  results 
provided  by  the  ode45  method. 


Figure  34.  Tension  for  a=TT/18  -  Q0= tt/10 
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Figure  35.  Maximum  and  minimum  tensions 


Figure  36.  Maximum  and  minimum  tensions  3-D  variation 
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Figure  37.  Minimum  tension  variation 
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Figure  38.  Maximum  tension  variation 
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In  Figure  34,  we  should  note  that  this  is  not  the  variation  of  tension  in  one 
of  the  strings,  but,  rather,  the  sum  of  them.  Although,  each  string’s  tension  is  zero 
for  every  other  half  period,  the  fact  that  we  describe  the  motion  of  the  mass  as  an 
equivalent  standard  pendulum,  with  altered  restoring  force,  and  we  do  not 
distinguish  the  angular  displacements  with  respect  to  the  different  strings  from 
each  other,  leads  us  to  treat  the  string  tension  with  the  same  methodology.  As  a 
result,  in  each  peak  in  Figure  34,  the  mass  switches  strings,  and  these  exact 
switches  take  place  when  the  string  tension  takes  its  maximum  value. 

From  Figures  35-38,  we  can  evaluate  the  variation  of  maximum  and 
minimum  tensions  exerted  on  the  strings  with  respect  to  both  a  and  0o.  On  the 
one  hand,  as  expected,  increasing  9o  increases  the  Tmax  and  decreases  Tmin. 
On  the  other  hand,  increasing  a,  and  keeping  6o  constant,  steadily  decreases 
Tmin.  An  interesting  point  is  that  Tmax  does  not  exhibit  the  same  behavior.  In 
Figure  39,  we  can  see  for  which  value  of  a  we  will  have  the  maximum  tension  for 
a  range  of  initial  displacements. 
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Figure  39.  Characteristic  angle  producing  the  maximum  tension  for  given  90 
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C.  NUMERICAL  SIMULATION 

A  very  precise  method  to  simulate  the  response  of  the  2-pendulum  is  with 
the  use  of  the  VISUAL  NASTRAN  4-D  software.  In  this  section  we  attempt  to 
validate  the  methods  that  we  used  to  predict  the  2-pendulum  response  in 
Chapter  II,  Section  A,  by  comparing  their  results  with  those  produced  by  VN  4-D. 

1 .  VISUAL  NASTRAN  Software  -  SOLIDWORKS 

To  simulate  the  response  of  the  2-string  pendulum  in  real  time,  a  specific 
procedure  had  to  be  followed.  This  procedure  was  initiated  by  designing  a  model 
of  a  2-string  pendulum  using  SOLIDWORKS.  This  3-D  design  software  allows 
the  exact  geometrical  design  of  structural  models,  which  are  compatible  with  the 
VN  4-D  software.  This  model  is  then  transferred  into  VN  4-D.  The  distinct 
attribute  of  VN  4-D  is  to  solve  differential  equations  in  3-D  space  and  time.  The 
final  model  used  in  VN  4-D  was  comprised  of  geometrical  parts:  mass  properties 
-  developed  in  SOLIDWORKS  — ,  and  their  assembly  -  developed  partially  in 
SOLIDWORKS  and  partially  in  VN  4-D.  The  assembly  constraints,  which  resulted 
in  the  rigidity  of  the  platform  where  the  2-pendulum  support  points  were 
positioned,  were  developed  in  SOLIDWORKS.  The  string  constraints  were 
developed  in  VN  4-D.  Before  initiating  each  VN  4-D  run,  the  mass  had  to  be 
positioned  in  the  exact  coordinates.  This  was  the  equivalent  to  a  combination  of 

initial  angular  displacement  O0  and  initial  angular  velocity  O0  of  the  VN  4-D.  The 

results  of  the  simulations  (in  terms  of  mass  position  wrt  which  is  the  model’s 
global  system  of  axis,  linear  velocity,  and  tension  of  the  strings)  were  then 
exported  into  Microsoft  Excel  spreadsheets.  In  addition,  as  defined  in  Figures  5 
and  6,  the  mass  position  results  were  processed  to  obtain  the  angular 
displacement  of  the  pendulum.  Finally,  all  the  above  results  were  graphed  with 
the  aid  of  Matlab.  This  allowed  comparing  them  with  our  analysis  results. 
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Figure  40.  Front  view  of  2-pendulum  VN  4-D  model 

Y 


Figure  41 .  Side  and  bottom  view  of  2-pendulum  VN  4-D  model 

The  front,  side,  and  bottom  views  of  the  model  that  were  used  for  the  VN 
4-D  simulations  can  be  visualized  in  Figures  40  and  41.  In  this  model,  one  can 


41 


adjust  the  mass  of  the  mass  (blue  sphere)  and  the  characteristic  angle  by 
altering  the  length  of  the  strings  (black  lines)  or  the  distance  between  the  support 
points  (red  cubes).  In  addition,  an  initial  displacement  can  be  applied  to  the  mass 
by  precisely  positioning  the  mass  to  a  specific  position  in  terms  of  an  exact  set  of 
coordinates  (X,Y  and  Z),  with  respect  to  the  Cartesian  global  system  of  axis 
represented  by  the  black  arrows.  It  should  be  noted  that  the  exact  positioning  of 
the  mass  is  not  possible.  For  a  few  of  the  first  time  steps  of  the  simulation,  this 
results  in  a  slight  deviation  of  the  expected  results.  The  numerical  integration 
parameters  that  were  used  in  all  the  VN  4-D  simulations  are: 


A/A 

Simulation  Parameter 

Value 

1 

Integration  Method 

Kutta-Menson 

2 

Animation  Time  Frame 

0.001  sec 

3 

Integration  Step 

Variable  -  0.0005  sec 

4 

Steps  per  Frame 

2 

Table  1.  VN  4-D  simulation  parameters 


The  above  parameter  combination  yielded  results  that  were  tabulated  in  a 
time-step  size  of  0.001  sec.  This  way  compatibility  was  achieved  with  the 
previously  presented  results  using  the  developed  Matlab  code. 

In  the  following  example,  we  chose  to  solve  the  case  of  the  2-string 
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pendulum  of  characteristic  angle  a=  —  and  initial  angular  displacement 


& 


— .  For  consistency  we  used  ^  =  lm,  while  the  actual  string  length  was 


/  =  1.0154m . 


The  VN  4-D  results  will  be  compared  to  the  predictions  of  the  three 
methods  described  in  Chapter  II,  Section  A,  numbers  6  and  7,  for  the  same 
values  of  t ,  l  ,a  and  60. 
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Figure  42.  VN  4-D  -  a=Tr/18  -  0o=tt/1O 
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Figure  43.  Theoretical  predictions  for  cr=Tr/18  and  @o=tt/10 
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Figure  45.  Typical  velocity  variation  for  2-pendulum 


2.  Evaluation  of  Results 

In  Figures  42  and  43,  we  can  see  the  oscillation  prediction  using  VN  4-D 
and  the  analytical  methods,  respectively.  To  more  precisely  appreciate  the 
discrepancies  between  those  methods’  predictions,  Figure  44  will  be  used.  We 
can  clearly  distinguish  from  Figure  44  that  VN  4-D  predicts  that,  with  time,  the 
amplitude  of  oscillation  decreases  and  the  angular  frequency  increases.  Further, 
for  angular  frequency,  the  comparison  was  indispensable.  We  decided  to  plot  the 
angular  frequency  every  half-cycle  to  be  able  to  capture  the  decrease  in  c o  in  the 
most  possible  detail.  From  Figure  44,  we  can  realize  that  ode45  provided  the 
most  accurate  prediction  for  the  go  for  the  first  half-cycle,  compared  to  SCPA  and 
Energy  methods.  Nevertheless,  after  that  point,  the  go  of  VN  4-D  started  to 
increase  almost  linearly.  This  increase  in  go  was  not  captured  by  the  three  other 
methods. 
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To  understand  why  this  increase  in  co  takes  place,  we  plotted  the  angular 
frequency  of  oscillation  predicted  by  VN  4-D.  As  seen  in  Figure  45,  for  every  half¬ 
cycle  when  the  mass  reached  the  point  where  it  switched  strings  about  which  it 
oscillated,  the  velocity  magnitude  dropped  by  a  certain  amount.  This 
discontinuity,  just  before  and  after  the  switch  point  in  velocity  magnitude,  raises 
questions  about  energy  conservation.  In  the  following  sections,  we  are  going  to 
investigate  in  more  details  the  dynamics  that  surface  at  this  point  and  their  effect 
in  terms  of  energy  loss  of  the  oscillation  mass. 

D.  THE  REAL  2-PENDULUM 

1.  Velocity  and  Force  Vector  Analysis-Impact  Dynamics 


Figure  46.  Velocity  and  force  vectors  for  2-pendulum  at  switch  point 

In  Figure  46,  we  can  visualize  the  forces  that  act  on  the  mass  during  this 
change  and  contribute  to  the  switch  of  center  of  rotation,  as  already  depicted  in 
Chapter  II,  Section  A,  number  1  and  Figure  7.  Nevertheless,  in  this  section  we 


46 


are  going  to  emphasize  the  way  that  the  forces  change  the  direction  of  the  linear 

velocity  from  Vb  to  Va,  as  well  as  other  effects,  such  as  energy  loss  and 
termination  of  oscillation  for  ar>45°. 

According  to  Figure  46,  the  dynamics  of  our  problem  strongly  resemble 
the  dynamics  of  a  ball  colliding  with  a  rigid  wall.  For  cr<450,  the  only  known 

parameters  are:  i)  vector  Vb  (direction,  orientation,  and  magnitude),  ii)  the 

direction  and  orientation  of  vector  Va,  iii)  the  weight  of  the  mass,  and,  finally,  iv) 
the  direction  and  orientation  of  7"/nsf.  Based  on  these  parameters,  we  will  attempt 
to  predict  the  response  of  the  real  2-string  pendulum. 

2.  The  Real  2-Pendulum  with  cr<45 

In  our  initial  analysis,  we  are  going  to  assume  that  the  string  is  stiff  enough 
and  allows  no  elongations  even  under  large  values  of  tension  forces.  Therefore, 
we  are  not  considering  displacements  along  the  direction  of  the  string. 
Nevertheless,  we  are  going  to  assume  that  an  amount  of  energy  is  lost  during 
this  switch  of  velocity  direction.  Following  the  discussion  of  Chapter  II,  Section  A, 
for  this  case  the  pendulum  will  oscillate  ideally  in  a  manner  that  was  already 
analyzed  in  that  section.  In  addition  to  the  previous  analysis,  we  are  going  to  take 
into  account  the  energy  losses  that  occur  at  the  switch  point  due  to  the 
introduction  of  nonlinearity.  For  this  purpose,  we  will  apply  the  equation  of  motion 
of  oscillation  between  two  arbitrary  consecutive  switch  point  impacts,  i  and  i  +  l , 

which  is  equation  (2.18),  0  +  ySin6,cos£Z  +  yCOs6,sinasign6,  =  O,  and  the 

simple  impact  rule  at  the  switch  point  i ,  for  instance.  The  simple  impact  rule 
states  that  the  impacts  are  instantaneous  and  the  energy  loss  would  be 
represented  by  a  reduction  in  the  velocity  magnitude.  That  is, 

etl=r0_,  (2.48) 

where  r  is  the  coefficient  of  restitution  at  impact  [3], 
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Referring  to  Figure  46,  we  could  be  more  specific  and  precise  while 
calculating  the  energy  absorption  during  impacts  described  by  equation  (2.48). 
As  already  mentioned,  the  impact  exerted  on  the  mass  by  the  string  can  be 
described  by  a  collision  of  a  rigid  ball  with  a  rigid  wall.  Since  the  weight  of  the 
mass  is  a  conservative  force,  it  does  not  absorb  any  energy  during  the  collision. 
Therefore,  Tjnst  is  the  only  force  that  contributes  to  energy  absorption,  which 
equals  the  work  done  by  force  W  during  the  period  of  contact  [4],  Moreover,  the 
sole  velocity  components  that  are  affected  by  W  are  those  parallel  to  the 
direction  of  Tinst.  Dealing  this  time  with  linear  velocities  (instead  of  angular),  we 

can  see  that  the  components  of  Vb  and  Va ,  which  are  parallel  to  the  direction  of 
the  string  causing  the  impact,  (heavy  dotted  line  in  Figure  46),  are  Vbag ,  with  a 


magnitude  of 


Vb 


aff 


Vb 


cos(90-2 a),  and 


Va 


aff 


Va 


cos(90-2  a).  Finally,  using 


Vb  =  V-i ,: Vaaff  =  Vaff+i  and  Vbaff  =  Vaff+i 


Vaaff,  with  a  magnitude  of 
the  notation  Va  =  V+j , 


JJ/  _  ^inst^Waff-j  +  Vgff +/)]  (2  49) 

where 

A Tm  =  Te_,  -  Tc+I  =  ^  =  pe_,  - pc_,  (2.50) 

and  At  is  the  infinitesimally  small  time  duration  of  the  impact. 

In  equation  (2.50),  Te_t  is  the  impact  force  exerted  by  the  mass  m  on  the 

string  causing  instant  elongation.  Tc+i  is  the  impact  force  exerted  by  the  string 
contracting  and  releasing  its  strain  energy  on  the  mass  m.  ,  pc_i  are  the 
respective  impulses.  Consequently,  in  the  circumstance  where  Tc+i  =  Te_f ,  there 

would  be  no  energy  loss  during  impact  and  we  would  deal  with  a  completely 
elastic  collision. 
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In  addition  to  the  above,  from  the  point  of  view  of  the  conservation  of 
momentum  and  considering,  once  more,  only  the  velocity  components  that  are 
parallel  to  the  string  direction,  during  the  phase  of  the  instantaneous  elongation 
of  the  string 


MV0  =  MVaff_l+Pel  (2.51) 

where  Vo  is  the  velocity  of  the  mass  at  the  end  of  the  elongation  phase,  pel  is 

the  reaction  impulse  force  that  is  exerted  on  the  mass,  and  M  is  the  effective 

mass  which  equals  M= - -.  In  this  case,  m  is  the  fictitious  mass  of  the 

m  +  m 

object  onto  which  the  mass  m  collides.  Since  this  object  is  represented  by  a  rigid 
wall  whose  momentum  is  infinite,  we  assume  that  m=  oo.  Therefore,  M  =  m 
and,  since  V0  -  0 ,  the  magnitude  of  the  impulse  force  is 

Pel=mVaff-i 

or 

pel  =  mV_t  cos(90-2cr)  =  Te_tAt  (2.52) 

where  At  is  the  infinitesimally  small  time  duration  of  the  impact  [3], 

The  next  step  should  be  to  calculate  the  pc  =  Tc  +  iAt.  If  we  conduct  a 
similar  analysis,  we  would  see  result 

pc  =  mV+i  cos(90-2a)  =  Tc+iAt  (2.53) 

Looking  more  closely  at  equation  (2.53),  we  realize  that  if  Tc+i  ^0  after 

the  string-change  point,  there  would  be  a  linear  velocity  component  parallel  to  the 
string’s  direction  in  addition  to  the  linear  velocity  component  perpendicular  to  the 
new  string’s  direction.  The  second  component  promotes  the  smooth  circular 
motion  of  the  mass,  while  the  first  component  would  push  the  mass  to  deviate 
from  the  circular  orbit.  Consequently,  if  we  want  to  predict  the  path  of  the  mass 
that  would  not  deviate  from  the  regular  pendulum’s  oscillation  path  we  will 
assume  that 
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(2.54) 


V+i  cos(90  -  2a)  =  0 

Thus,  equation  (38)  should  be  manipulated  and  corrected  as  follows 

p  _V+i_  l(Ki  sin(90  -  2a))2  +  (V+i  cos(90  -  2a))2  _ 

+‘  l  \  r- 


( V_j  sin(90  -  2a))2  +  (rV_t  cos(90  -  2a))2 

\  T-  ” 

TO sin(90  -  2a))2  +  (r  cos(90  -  2cr))2  <^> 

=  6L;  ^/(sin(90  -  2cr))2  +  (r  cos(90  -  2cr))2  (2.55) 

where,  considering  (2.55),  the  coefficient  of  restitution  should  be  r  =  0  . 

The  cancellation  of  the  parallel  component  of  the  velocity  during  impact  is 
not  only  an  assumption  that  facilitates  the  work  of  predicting  the  response  of  a 
real  2-pendulum,  but,  also,  a  very  good  approximation.  The  sudden  impulse 
shock  that  is  applied  on  the  string  from  the  mass  produces  a  compressive  wave 
that  propagates  throughout  the  length  of  the  string.  This  wave  will  finally  die  out 
under  the  effect  of  damping.  Damping  in  almost  inextensible  strings  may  be 
caused  by  several  factors,  such  as  viscosity  of  the  air  and  internal  friction.  The 
velocity  of  elastic  wave  propagation  within  the  string  mass  is  similar  to  the  speed 
of  sound  [5] [6].  Therefore,  by  comparing  this  velocity  with  the  time  scale 
associated  with  the  periodic  motion  of  the  pendulum,  the  time  required  for  the 
completion  of  the  energy  loss  related  to  this  shock  is  negligible.  All  the  above 
support  the  validity  of  the  assumption  that  the  coefficient  of  restitution,  associated 
with  the  impact  of  the  mass  on  the  string,  is  zero. 

Finally,  from  equations  (2.54)  and  (2.55),  the  amount  of  energy  loss  at  the 
switch  point  is 

=^mV2^  ((1  -  r)  cos(90  -  2a))2  (2.56) 
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The  system  of  equations  (2.18)  and  (2.48)  or  (2.55)  can  be  solved 
numerically  using  a  constant  time  step  of  A?  =  0.001sec  and  the  forward 
difference  scheme.  This  takes  into  account  equation  (2.48)  or  (2.55).  The  Matlab 
code,  Program_3.m,  provides  the  solution  to  the  system  of  equations  for 
predetermined  number  of  periods. 

In  the  following  sections,  we  are  going  to  compare  the  prediction  of  the  2- 
pendulum  response  in  terms  of  angular  displacement,  angular  velocity,  and  string 
tension.  We  will  use  the  above  m-file  (solving  equations  (2.18)  and  (2.55)) 
combined  with  those  obtained  from  VISUAL  NASTRAN.  Additionally,  we  will  plot 
the  phase  diagrams  for  each  method,  as  well  as  the  discrepancy  between  the 
prediction  and  the  VISUAL  NASTRAN  simulation  results.  This  will  better  validate 
the  precision  of  the  Matlab  code. 
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Figure  47.  VISUAL  NASTRAN  -  cr= tt/1  8  -  0o=Tr/1 0 
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Figure  48.  MATLAB  -  gftt/1  8  -  0o=tt/1  0 
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Figure  49.  Prediction  discrepancy 
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Figure  50. 
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Figure  51.  MATLAB  -  a=Tr/18  -  0o=Tr/1O 
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Figure  52. 
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Figure  53.  VN  4-D  -  Phase  diagram  -  cr=Tr/18  -  0o=tt/1O 
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Figure  55.  VISUAL  NASTRAN  -  cftt/1  2  -  Q0=t r/1 0 
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Figure  57.  Prediction  discrepancy 
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Figure  58.  VISUAL  NASTRAN  -  a=ul  1 2  -  0o=tt/1  0 
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Figure  59.  MATLAB  simulation  -  a=TT/12  -  Q0= tt/10 
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Figure  60. 


Prediction  discrepancy 
58 


9'(rad/sec)  0'  (rad/sec) 


1.5 


-0.3  -0.2  -0.1  0  0.1  0.2  0.3 


6  (rad) 

Figure  61 .  VISUAL  NASTRAN  Phase  diagram 
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Figure  62. 
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Figure  63.  VISUAL  NASTRAN  -  a=TT/6  -  0o=tt/1  0 
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Figure  64.  MATLAB  simulation  -  a=Tr/6  -  @0=tt/10 
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Figure  66. 
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Figure  67.  MATLAB  simulation  -  ar=Tr/6  -  0o=Tr/1O 


t  (sec) 

Figure  68.  Prediction  discrepancy 

62 


0' (rad/sec)  9'  (rad/sec) 


2 


1.5 

1 

0.5 

0 


-0.5 


-1.5 


1  1 

1  1 

1  1 

1  1 

'll  1 

1  1 

1  1 

i _ f 

X  ! 

L\ . L. 

T  i 

/  ; 

/  i 

i 

i 

i 

i 

i 

i 

i  i  / 

;  ^ 

1  i 

-2 


--I - I - L  - 

-0.3  -0.2  -0.1 


■i - 1 - 1 - 1-  . 

0  0.1  0.2  0.3 


9  (rad) 


Figure  69.  Phase  diagram 


Figure  70.  Phase  diagram 
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Evaluating  the  above  set  of  figures,  we  conclude  the  following: 

•  The  variable  frequency  of  oscillation,  as  well  as  when 
the  string  switch  takes  place,  is  relatively  well  predicted.  The  error  between  the 
predicted  and  actual  performance  causes  spikes  in  the  angular  velocity 
discrepancy  plots. 

•  In  terms  of  the  prediction  of  the  angular  displacement 
and  velocity  of  the  oscillation,  the  discrepancies  are  caused,  on  the  one  hand,  by 
the  failure  to  precisely  predict  both  the  location  in  time  of  the  switch  points  and 
the  energy  loss  at  every  one  of  them.  On  the  other  hand,  it  was  not  possible  to 
obtain  the  exact  equations  used  in  VN  4-D  to  calculate  these  parameters. 
Therefore,  a  comparison  with  the  VN  4-D  software  results  does  not  provide  a 
precise  appreciation  of  the  actual  accuracy  of  the  results. 

•  For  smaller  values  of  characteristic  angle  a,  which 
equals  less  energy  loss  at  the  switch  point,  the  prediction  discrepancy  is 
minimized. 

b.  String  Tension  Comparison 

For  the  string  tension  prediction,  equation  (2.46)  was  used.  The 
Matlab  code  had  already  computed  the  angular  velocity  and  displacement  inputs. 
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Figure  7 1 .  Tension  comparison  a=nY1 8  -  Q0= tt/1  0 


Figure  72.  Tension  comparison  a=Tr/12  -  @o=Tr/10 
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Figure  73.  Tension  comparison  a=Tr/6  -  @o=tt/10 

The  above  plots  lead  to  the  following  conclusions: 

•  It  can  be  seen  that,  for  larger  values  of  a,  which  correspond 
to  extensive  energy  loss  at  the  switch  point,  VN  4-D  reveals  time  periods  when 
the  tension  on  both  strings  is  zeroed.  This  corresponds  to  an  instantaneous 
jump.  Because  the  mechanism  of  energy  restitution  of  the  ropes  used  in  VN  4-D 
is  not  known  and  that  there  is  the  possibility  that  the  previously  assumed  zero 
parallel  velocity  component  is  not  completely  eliminated,  the  resulting  jump  after 
impact  could  be  expected. 

•  The  tension  predicted  by  VN  4-D  shows  spikes  of  irregular 
magnitude  at  the  switch  points.  The  existence  of  these  features  can  be  related  to 
the  previously  presented  explanation.  Although  they  do  not  show  a  pattern,  the 
fact  that  a  previously  slack  string  becomes  taught  within  an  infinitesimally  small 
time  duration  can  provoke  this  peak  tension  magnitude.  This  corresponds  to  the 
impulse,  which  zeroes  the  parallel  velocity  component. 
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•  The  Matlab  code  tension  prediction  coincides  with  the  VN  4- 
D  only  at  the  instances  when  the  mass  velocity  is  zero  (indicated  with  a  purple 
circle  in  Figure  72).  Therefore,  the  VN  4-D  tension  prediction  is  not  in  accordance 
with  equation  (2.46).  To  understand  this  unexpected  result,  a  model  of  two 
standard  pendula  ~  one  with  a  string  and  the  other  with  a  massless  rod  --  was 
used  to  calculate  the  tension  during  free  oscillation  with  an  arbitrary  initial 
displacement.  The  results  revealed  that  the  standard  rod  pendulum’s  tension 
agrees  with  equation  (2.46),  while  the  standard  string  pendulum’s  tension  does 
not.  Moreover,  an  effort  was  made  to  understand  the  pattern  that  provided  the 
string  tension  results.  After  some  trial  and  error  efforts,  the  equation  that  provided 
an  almost  identical  tension  prediction  was 


rr  - m02l  ,  n 

T  =  — - — +  mg  cost) 


(2.57) 


Efforts  were  made  to  understand  the  above  equation,  but,  unfortunately, 
the  MSC  Software  no  longer  supported  the  specific  software.  Consequently, 
further  tension  comparisons  with  VN  4-D  will  not  be  performed  and,  also,  the 
Matlab  results  will  be  presented  without  validation  from  the  VN  4-D  results. 


3.  The  Real  2-Pendulum  with  a>45 

In  the  case  that  ar=45°,  when  applying  equation  (2.55),  we  realize  that  for 
r  =  0  all  the  kinetic  energy  of  the  pendulum  would  be  absorbed  by  the  string 
during  impact.  Thus,  the  oscillation  would  be  terminated.  The  above  can  be 
visualized  from  the  following  graphs,  where  VISUAL  NASTRAN  and  MATLAB 
simulations  are  used  to  provide  the  response  of  the  2-string  pendulum  with 
a=45°. 
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Figure  74.  VISUAL  NASTRAN  -  a=n/4  -  Q0= tt/1  0 
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Figure  75.  MATLAB  simulation  -  or=u/4  -  0o=Tr/1O 
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Figure  76.  Prediction  discrepancy 
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Figure  79.  Prediction  discrepancy 
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Figure  80.  VISUAL  NASTRAN  -  a=Tr/4  -  0o=Tr/1O 
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Figure  82.  Tensions  cr=Tr/4 

Although  the  limiting  case  of  a=45°  can  be  described  by  equations  (2.18) 
and  (2.55)  (assuming  only  that  r  =  0  )  if  a>45°,  the  motion  of  the  mass  will  not  be 
governed  by  linear  or  nonlinear  equations  of  motion,  but  only  by  its  interaction 
with  the  strings. 


Force  and  velocity  diagram  for  a>45° 
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Figure  83. 


According  to  Figure  83,  the  linear  velocity  before  the  switch  point  Vb  (red 

vector)  is  composed  of  two  components,  Vbi  and  Vbi  (turquoise  and  green 
vectors,  respectively).  At  the  instance  when  the  right  string  becomes  tight,  the 

tension  of  the  string  (blue  vector)  exerts  an  impulse  affecting  Vbi .  Irrespective  to 

the  value  of  the  restitution  coefficient  that  governs  this  impact,  Vbi  cannot 
promote  the  circular  motion  of  the  mass  with  respect  to  the  right  string  due  to  its 
direction.  Contrarily,  the  left  string,  which  becomes  slack  after  the  switch  point  for 
cr<45°,  continues  to  be  tight  and  exerts,  along  with  the  right  string,  an  impulse  on 

the  mass  opposing  Vbi  (orange  vectors).  Thus,  either  the  mass  will  bounce  or  it 
will  come  to  rest.  This  depends  on  the  restitution  coefficient.  Consequently,  if  we 
preserve  the  assumption  that  r  =  0  for  all  mass-string  impact  interactions,  after 
the  mass  reaches  the  switch  point,  it  comes  to  rest. 

E.  THE  DRIVEN  PENDULUM 

Extending  the  analysis  performed  in  the  previous  sections,  we  should 
include  the  case  of  the  pendulum  in  which  the  supporting  point,  or  points, 
undergoes  motion  due  to  forcing.  The  moving  support  induces  complexity  to  the 
response  of  the  pendulum  that  may  not  end  in  oscillatory  behavior  [6], 

1.  The  Standard  Driven  Pendulum 

Consider  a  standard  pendulum  which,  at  an  arbitrary  moment  of  the 

7Z  TC 

oscillation,  has  an  angular  displacement  6,  where  —  .  We  define  a 

world  system  of  axis  where  x-axis  points  to  the  right  and  y-axis  points  upwards. 
An  arbitrary  motion  drives  the  support  point  of  the  pendulum,  whose  acceleration 
at  that  moment  has  both  a  horizontal  and  a  vertical  component. 
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Figure  84.  The  standard  driven  pendulum 


According  to  the  previous  analysis,  the  equation  of  motion  for  the  undriven 
standard  pendulum  is  either  equation  (2.6)  or  (2.2).  The  effect  of  the  horizontal 
component  of  the  support  acceleration  is  an  inertial  force  on  the  mass 

Fx  =  -mX  .  Similarly,  for  the  vertical  component,  the  inertial  force  is  Fy  =  -mY  . 
We  notice  that  Fy  has  the  same  direction  and  orientation  as  gravity.  These  two 

forces  contribute  to  restoring  force  with  two  additional  components  depending  on 
angle  9.  Therefore,  the  new  restoring  force  is 

Fr  =  (mg  +  mY)sin0-mXcos0  (2.58) 

Therefore,  the  equation  of  motion  becomes 

0  +  ^  +  ^  sin  0  -  cos  =  0  (2.59) 
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Finally,  from  equation  (2.59)  we  can  draw  the  conclusion  that,  due  to  the 
horizontal  acceleration,  the  pendulum  will  initiate  oscillations  even  without  the 
existence  of  initial  displacement  or  initial  velocity.  Thus,  for  6=Q0= 0,  equation, 

X 

(2.59)  becomes  0  =  —  .  Finally,  the  above  equation  is  valid  for  the  range  of  9 

already  mentioned  with  one  exception:  when  a  jump  occurs.  The  jump  behavior 
of  the  driven  pendulum  will  be  analyzed  in  a  subsequent  section. 

2.  The  Driven  2-Pendulum 

Based  on  the  previous  analysis  and  equation  (14),  we  can  derive  the 

71  71 

equation  of  motion  for  the  driven  2-pendulum  as  follows  for  <0+a<— : 


6  +  sin(#  +  a)-^~  cos(0 +a)  =  0 


0 


(sin  $  cos  a  +  signO  cos  0  sin  ( signO  cos  0  cos  a  -  sin  0  sin  a)  =  0 


0  +  (^~y —  cos  a  +  -j-  sin  O')  sin  6  +  sin  a  -  cos  a)signO cos  6  =  0  (2.60) 

Again,  the  validity  of  equation  (2.60)  is  questioned  when  a  jump  occurs. 
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Figure  85.  The  driven  2-pendulum 

a.  String  Tension  Investigation  of  the  Driven  2-Pendulum 

Assuming  that  the  support  of  the  2-pendulum  is  accelerated  in  the 
direction  shown  in  Figure  (85),  the  mass  will  leave  the  equilibrium  position  only  if 
T2=0  N.  For  this  to  occur,  the  vertical  component  of  T1  should  be  able  to 

withstand  the  virtual  weight  of  the  mass,  which  is  Wv  =  m(Y  +  g)  .  We  can  easily 

prove  that  the  magnitude  of  the  string  tension  about  which  the  mass  will  oscillate 
is 

T  =  m(02l  +  X  sin(9 + a) +  Ycos(9  + a) +  gcos(0  + a))  (2.61) 

7T  7Z 

which  is  valid  for  all  +  a<G<—  -a .  Nevertheless,  close  attention  should 

2  2 

be  given  in  terms  of  the  sign  of  angles  Q  and  a.  If  inconsistencies  are  to  be 

avoided,  angle  o  should  always  have  the  same  sign  as  angle  6.  A  quick  look  at 
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Figure  86  shows  us  that  a  positive  X  would  increase  the  tension  on  the  string 
while  a  negative  X  would  decrease  it.  Nevertheless,  if  the  mass  were  oscillating 

about  the  other  string,  a  negative  X  would  increase  the  string  tension,  and  so 
on.  The  above  facts  are  predicted  by  equation  (2.61 ). 


Figure  86.  2-pendulum’s  force  and  acceleration  diagram  during  oscillation 

Nevertheless,  equation  (2.61)  has  one  limitation:  the  tension  can 
only  take  positive  values.  When  the  string  tension  is  zeroed,  a  jump  occurs.  That 
behavior  will  be  analyzed  in  a  subsequent  section. 

In  addition,  it  would  be  useful  to  know  the  magnitude  of  the  tension 
exerted  on  the  strings  while  the  mass  remains  in  its  equilibrium  position.  While 
the  mass  is  in  that  position,  the  sum  of  the  vertical  components  of  T1  and  T2 
should  be  equal  to  the  virtual  weight  of  the  mass  Wv,  while  the  difference  of  their 

horizontal  components  should  be  equal  to  mX 
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cos a(Tl  +  T2)  =  m(Y+g)  and  sina(Tl-T2)  =  mX 

solving  the  above  system  of  two  equations 

n_mAY  +  g)  |  X  , 

2  cos  a  sin  a 

Y2  =  —  +  s)  _ 

2  v  cos  a  sin  a J 


(2.62) 


b.  Stability  of  the  2-Pendulum 

While  equation  (2.60)  describes  the  2-pendulum’s  oscillation  in  the 
absence  of  non-zero  initial  conditions,  the  oscillation  will  not  be  initiated  under 
any  value  of  horizontal  acceleration. 


Manipulating  equation  (2.60)  with  zero  initial  conditions,  results  in 
9  =  -  ^ y  ^  sin  a  +  ^-cos  a .  The  vertical  components  of  the  restoring  force 

tend  to  drive  the  mass  to  move  downwards.  This  motion  is  constrained  in  this 
situation  since  the  mass  is  supported  by  both  the  taught  string  (an  exception 

exists  for  the  case  that  Y  is  negative  and  g  +  F<0).  Therefore,  an  angular 


acceleration  will  be  induced  only  if  the  horizontal  component  manages  to 
overcome  the  condition  and 


is  fulfilled. 


X 

tana 


(2.63) 


Another  approach,  which  would  lead  us  to  the  same  conclusion,  is 
analyzing  the  string  tension.  Assuming  that  the  support  of  the  2-pendulum  is 
accelerated  in  the  direction  shown  in  Figure  86,  the  mass  will  leave  the 
equilibrium  position  only  if  T2=  0  N.  For  this  to  occur,  the  vertical  component  of  T1 
should  be  able  to  withstand  the  virtual  weight  of  the  mass,  which  is 
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Wv  =  m(Y+g).  In  the  case  of  the  2-pendulum,  which  is  about  to  leave  its 

equilibrium  position  0  =  0  and  from  (2.61),  for  the  mass  to  be  lifted  up  we  must 
have 

T  cos  a  =  m{X  sin  a  +  Y  cos  a  +  g  cos  a)  cos  a  >  m(Y  +  g) 

X  sin  a  cos  a  >  g  sin2  a  +  Y  sin2  a 


tan  a 

Finally,  from  equation  (2.63),  we  can  draw  the  conclusion  that 
choosing  larger  values  for  a  would  make  the  2-pendulum  more  stable. 
Nevertheless,  this  advantage  has  a  considerable  drawback:  the  increase  in  string 
tension  associated  with  larger  characteristic  angles.  In  Figure  87,  one  can 
observe  both  the  limiting  value  of  horizontal  acceleration,  which  allows  the  2- 
pendulum  equilibrium,  and  the  maximum  tension  associated  with  this  horizontal 
acceleration  for  a  range  of  characteristic  angles  (1°<a<80°  ). 
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Figure  87.  Horizontal  acceleration  inducing  instability  versus  associated  maximum 

tension 
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From  the  above  figure,  it  can  be  seen  that,  although  the  necessary 

X  would  break  the  equilibrium,  which  increases  steadily  for  the  first  half  of  the 
range  of  cr,  the  related  tension  increases  in  a  slower  pattern.  This  difference  in 
increasing  pattern  is  caused  by  the  analytical  equations  that  describe  them. 
These  equations  are 


Xlim=(S  +  7min)tan« 

T 


_  (g  +  y max)  ! 


(2.64) 


cos  a  sm  a 

One  method  to  evaluate  the  optimum  choice  of  a  would  be  to  define 
a  ratio  of  acceleration  over  the  maximum  attained  tension.  According  to  the 
previous  equations,  this  ratio  would  be  R  =  sin  a . 


Nevertheless,  this  method  is  not  effective.  The  values  of  the 

required  X  to  disturb  equilibrium  are  very  large  for  the  second  half  of  the  a 
range,  and  even  for  some  values  of  the  first  half.  A  more  realistic  approach  is  to 
choose  a  value  of  a  (for  example  20°)  as  a  reference  value,  aref,  and,  then, 
compare  the  maximum  tension,  T maxref ,  with  the  one  exerted  on  the  strings  of  the 

2-pendulums  of  larger  cr,  but  still  associated  with  the  same  Xrej- . 
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Figure  88.  a=20°  reference  angle 

The  above  plots  reveal  that,  in  a  choice  of  a  greater,  but  less  than 
52°,  a  would  be  a  better  choice. 

Another  approach  is  to  find  the  angle  a  which  will  result  in  the 
minimum  exerted  tension  for  a  combination  of  vertical  and  horizontal 
accelerations.  In  the  following  figures,  we  can  visualize  the  variation  of  optimum 

o  when  X  and  Y  vary  from  zero  to  g  along  with  associated  tension.  The  results 
are  validated  with  the  constraint  that  the  2-pendulum  should  be  stable  for  the 

combination  of  a,  X  ,  and  Y . 
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Tension  (N)  CD  a  (degrees) 


Characteristic  angle  giving  the  minimum  tension 


89.  Optimum  characteristic  angle  investigation  for  combinations  of 

accelerations 


Minimum  tension  associated  with  optimum  characteristic  angle 


Figure  90.  Tension  related  to  optimum  a 
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A  first  step  to  take  advantage  of  the  above  plots  and  analysis  would 

be  to  choose  the  maximum  expected  X  and  Y  for  our  design  and,  then,  iterate 
for  the  optimum  a  in  terms  of  minimum  tension.  The  second  step  would  be  to 

verify  that  the  design  would  be  stable  for  a  new  range  of  X  and  Y  composed  of 
negative  values  because  Y  tends  to  destabilize  the  2-pendulum.  In  the  following 
example,  the  range  of  -0.3g<7  <0.7 g  is  chosen  to  challenge  the  stability  of 
the  selection  of  a=38°  designed  to  provide  the  minimum  tension  for  the 
X  =  g,Y  =  g  combination  of  maximum  accelerations. 

As  depicted  in  the  following  figure,  there  exists  a  range  of  X  and 
Y ,  where  the  driving  design  constraint  is  not  the  desired  minimum  tension,  but, 
rather,  the  stability  for  the  pendulum.  Therefore,  for  X  =  Y  =  g  ,  we  should  select 
another  characteristic  angle  according  to  the  stability  criterion. 


Design  Driving  Constraint 


Figure  91 . 


Design  driving  factor  constraint  versus  objective  function  minimization 
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Characteristic  angle  for  combined  constraints 


Figure  92.  Characteristic  angle  selection  for  combined  constraint  and  objective 

function  minimization 


Therefore,  from  Figure  92,  we  can  decide  which  should  be  angle  a. 
This  would  be  based  on  the  worst  expected  negative  vertical  acceleration.  For 

example,  a  combination  of  X  -  g,Y  =  -0.3 g  would  require  a= 55°. 


In  addition  to  the  above  manipulations,  another  design  method 
would  be  the  use  of  basic  optimization  techniques.  In  terms  of  optimization,  the 
problem  is  a  minimization  problem,  i.e.,  minimizing  the  maximum  exerted  tension 

(71ax  -  ^  +  ^liax^  +  ^max  )  with  one  objective  variable  (characteristic  angle  o) 
cos  a  sincf 


X 

and  subject  to  boundary  conditions  (a  >  0,a  >  arctan( - ^ — )  ).  A  Matlab  code 

g  +  Y  ■ 

o  min 

was  developed  to  yield  the  optimum  characteristic  angle  for  given  maximum  and 
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minimum  accelerations.  In  the  following  figure,  we  can  visualize  the  results  for  an 
arbitrary  input  combination. 


maxXddot=0.4g  -  maxYddot=0.3g  -  minYddot=-0.3g 


Figure  93.  Optimum  selection  for  arbitrary  design  inputs 

c.  Jumps  of  the  Driven  2-Pendulum 

The  tension  magnitude  provided  by  equation  (2.61)  denotes  that, 
under  certain  conditions,  string  tension  may  become  zero.  Keeping  the 

71  TC 

assumption  that  -—  +  a<6<—-a  ,  the  only  terms  that  can  take  tension  off  the 

string  are  Xsin#  and  Y cos#.  Therefore,  the  restriction  that  ||#  +  a||>—  for  a 

jump  to  occur,  which  is  valid  for  the  unforced  2-pendulum,  is  unnecessary.  To 
predict  the  motion  of  the  mass  after  the  jump,  we  should  follow  a  certain 
procedure. 
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Jump 

Jump 

end 

initiation 

Figure  94.  Jump  of  the  2-pendulum 


Firstly,  just  after  the  jump,  we  should  track  the  position  with  respect 
to  a  global  coordinate  system,  x0  and  y0,  and  to  the  linear  velocity  magnitude,  u, 
of  the  mass.  With  the  above  initial  conditions,  the  mass  will  follow  a  projectile 
path.  Thus,  the  equations  of  motion  are 

x  =  x0  +  u0t  cos  </) 

'  1  2  (2.Q5) 

y  =  y0  +  u0tsm(/>--gt- 

Secondly,  after  the  jump,  we  should  track  the  positions,  velocities, 
and  accelerations  of  the  two  support  points.  Depending  on  the  type  of  motion 
applied  to  the  support,  we  should  especially  track  the  position  of  the  support 
points  during  the  jump  for  each  time  step. 

Finally,  we  should  track  the  distance  between  the  mass  and  each 
one  of  the  support  points.  When  this  distance  becomes  equal  to  the  string  length, 
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it  is  dynamically  possible  for  the  oscillation  to  proceed  again.  The  methodology 
used  to  predict  the  oscillation  characteristics  after  the  end  of  the  jump  is  similar  to 
the  one  used  while  analyzing  the  pendulum’s  behavior  at  the  switch  points. 
Firstly,  the  jump  takes  place  while  both  the  strings  lengths  are  less  than  / . 
Consequently,  the  first  step  is  to  identify  which  string  becomes  taught  first.  This 
determines  if  angular  displacement  has  occurred  after  the  jump.  Secondly,  the 
magnitude  of  post-jump  angular  displacement  has  to  be  calculated.  Knowing  the 
mass  positions  just  before  the  end  of  the  jump,  x  and  y ,  as  well  as  each  of  the 

string’s  support  point  positions,  Zleft  =  Z-/sin  a  and  Zri  ht  =  Z  +  /  sin  a  ,  the 


total  angle  defined  by  the  taught  string  and  the  vertical,  cp  is 


(/)  =  tan(  \"ght  =.),  0>O 


(f)  =  tan( 


z-Z 


(2.66) 


left 


yfieft—y 


),0<O 


Therefore, 


0  =  (j)-a,6  >  0 
6  =  -(f)+a,6  <  0 


(2.67) 


Moreover,  at  the  moment  when  one  of  the  strings  becomes  taught 
one  faces  again  a  mass  to  string  impact.  Knowing  the  linear  velocities  of  the 
mass  just  before  the  impact, 


ux  =  u0  cos  (j) 

Uy  —  Uq  sin  (j)-gt 


(2.68) 


and  the  horizontal  -  vertical  velocities  of  the  support,  X  and  Y ,  the  mass  will 
undergo  impact  with  the  taught  string.  This  is  governed  by  the  same  restitution 
coefficient  assumption  that  was  explained  in  Chapter  II,  Section  D,  number  2.  In 
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this  case,  the  impact  velocity  components  are  the  relative  velocities  between  the 
mass  and  the  support,  Uximp  =  ux - X,Uyimp  =uy-Y .  Consequently,  the  following 
two  equations  are  used  to  obtain  the  post-jump  angular  velocity 


9  =  ) {~Uximp  cos  (/)  +  Uyimp  sin  </>),  6  >  0 
9  =  y  (~Uximp  cos ^ - Uyimp  sin  (/)),0<  0 


(2.69) 


Finally,  knowing  the  post-jump  angular  displacement  and  velocity 
initial  allows  the  continuation  of  the  oscillation  subject  to  the  new  initial 
conditions. 


3.  Evaluation  of  Results 

In  the  following  two  sets  of  figures  (Figures  95-99  and  Figures  100-104), 
two  examples  of  the  response  prediction  for  forced  oscillation  with  jumps  are 
presented.  A  Matlab  code  Program_3.m  was  developed  to  predict  the  position  of 
the  mass  while  the  supports  are  subjected  to  horizontal,  for  the  first  set,  and  to 
horizontal  and  vertical  motion,  for  the  second  set,  under  conditions  that  will  force 
a  jump  to  occur.  The  motion  that  was  applied  to  the  support  midpoint  can  be 
visualized  in  Figures  95  and  100.  It  can  be  seen  that  in  order  to  produce  a  jump, 
i.e.,  zero  tension  on  both  strings,  we  applied  a  sudden  stop  to  a  motion  that  was 
previously  under  the  effect  of  variable  acceleration.  This  violent  stop  was  enough 
to  produce  a  distinguishable  jump  both  in  VN  4-D  and  in  Matlab.  The  blue  solid 
lines  represent  the  Matlab  predictions  for  the  mass  positions  and  angular 
displacement,  while  the  red  dashed  lines  represent  the  VN  4-D  predictions. 
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Figure  95.  Support  midpoint  trajectory 
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Figure  96.  Mass  trajectory 
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absolute  vertical  position  (m)  absolute  horizontal  position  (m)  0(rad) 


Figure  97.  Angular  displacement 


Figure  98.  Mass  absolute  position 
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Figure  99. 


Tension  variation 


2  i  i  i  i  i  i  i  i  i  i 

1.5 

1 

0.5 

0 

0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5  5.5 

t  (sec) 

2 
1.5 
1 

0.5 
0 

0  0.5  1  1.5  2  2.5  3  3.5  4  4.5  5  5.5 

t  (sec) 

Figure  100.  Support  midpoint  trajectory 
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global  horizontal  position 


t  (sec) 


Figure  101 .  Mass  trajectory 


Figure  102.  Angular  displacement 
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absolute  vertical  position  (m)  absolute  horizontal  position  (m) 
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Figure  103.  Mass  absolute  position 


Figure  104.  Tension  variation 
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On  the  one  hand,  the  positions  of  the  support  midpoint  (Figures  95 
and  100)  and  the  total  mass  positions  (Figure  96  and  101)  are  calculated  to 
coordinate  the  axis  system  with  its  origin  positioned  on  the  support  midpoint  Xg- 
Yg.  On  the  other  hand,  the  absolute  mass  position  is  calculated  wrt  a  coordinate 
axis  system  with  its  origin  at  the  equilibrium  position  of  the  mass  following  the 
support  motion  (Figures  98  and  103).  Moreover,  the  angular  displacement  is 
presented  where,  during  the  jumps,  the  value  is  set  to  zero  (Figures  97  and  102). 
Finally,  as  explained  in  Chapter  II,  Section  D,  number  2.d,  the  strings’  tensions 
are  plotted  without  including  the  VN  4-D  results. 

The  discrepancy  that  can  be  observed  between  the  Matlab  and  VN 
4-D  results  does  not  seem  to  be  related  to  the  complexity  of  the  induced  support 
motion.  After  a  more  extensive  study  of  the  VN  4-D  results  (tabulated  data 
exported  from  the  software  interface),  it  can  be  seen  that,  in  some  of  the  jumps,  a 
nonlinear  change  in  the  mass  velocity  exists  at  the  exact  moment  of  the  jump 
initiation.  Introducing  this  exact  nonlinearity  in  the  Matlab  code  yields  the  same 
results  with  VN  4-D. 

4.  Convergence  Evaluation 

An  essential  and  final  step  that  must  be  performed  before  ending  the 
analysis  of  the  driven  2-pendulum  is  to  examine  the  convergence  of  the  results 
for  different  choices  of  time  steps.  From  the  initial  analysis  of  the  free  2- 
pendulum  to  the  simulation  of  the  driven  2-pendulum  with  jumps,  the  integration 
time  step  selected  for  the  numerical  solution  of  the  nonlinear  equations  of  motion 
was  c/f=0.001  sec.  Therefore,  a  basic  evaluation  procedure  of  the  examples  will 
be  presented  to  verify  that  the  selected  time  step  offers  convergence  of  results 
over  other  choices  of  larger  time  steps. 
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Figure  1 05.  Convergence  evaluation  for  forced  oscillation  -  no  jumps 
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Figure  106. 


Convergence  evaluation  for  forced  oscillation  -  jumps 
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In  Figure  105,  the  convergence  is  tested  with  both  vertical  and  horizontal 
support  motion.  A  case  with  jumps  is  presented  in  Figure  106.  From  both  figures, 
one  can  easily  conclude  that  gradually  increasing  the  time  step  dt  will  provoke  a 
solution  divergence.  While  the  discrepancy  for  dt= 0.005  sec  or  dt=  0.01  sec  (up  to 
ten  times  larger  than  the  initial  choice)  will  cause  the  solution  to  change  slightly,  a 
further  increase  in  the  time  step  size  may  yield  solutions  with  larger 
discrepancies  than  the  original.  Especially  in  the  jump  cases,  it  can  be  seen  that 
some  jumps  are  incorrectly  predicted  (in  terms  of  time  of  jump  initiation  and 
termination)  or  even  skipped  if  the  time  step  is  large  enough. 
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III.  THE  FOUR-STRING  PENDULUM 


A.  THE  4-PENDULUM  DESIGN  ANALYSIS 


To  extend  the  previous  analysis  of  the  2-pendulum  to  more  complex 
designs,  the  performance  of  a  four-string  pendulum  will  be  analyzed.  The  four 
strings  have  the  same  length  and  their  supports  are  symmetrically  placed  wrt  the 
mass.  All  the  strings  are  attached  to  the  center  of  mass.  The  following  figures 
represent  the  model  that  was  used  in  VN  4-D  and  offers  four  views  of  the  design: 


Figure  108.  4-pendulum  bottom  and  side  view 
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Such  a  design  is  one  of  the  more  representative  simplifications  of  the  real 
world  suspension  designs.  Multi-rope  suspension  designs  differ  from  each  other 
on  the  number  of  ropes  that  are  used,  the  positioning  of  their  support  points,  and 
the  selection  of  their  attachment  positions  on  the  suspended  load. 

The  purpose  of  this  section’s  analysis  will  not  be  to  investigate  the 
swinging  attitude  of  such  a  pendulum,  but,  rather,  to  study  its  ability  not  to  swing 
when  accelerating  motion  is  applied  to  its  support.  Moreover,  the  string  tensions 
will  be  analyzed  to  try  to  draw  the  most  accurate  picture  of  the  design 
characteristics  that  should  be  considered  to  achieve  a  desirable  4-pendulum 
design. 


Figure  109.  4-pendulum  design  characteristics 


In  the  above  figure,  the  black  sphere  represents  a  Ikgr  mass  and  the 

black  lines  are  the  strings  of  the  pendulum.  The  global  coordinate  system  of  axis 

is  presented.  It  dictates  that  the  positive  X-axis  represents  the  forward  motion; 
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the  positive  Z-axis  represents  the  side  motion  towards  the  right;  and  the  positive 
Y-axis  represents  the  upward  motion  of  the  support.  Since  we  are  now  dealing 
with  a  3-D  analysis,  instead  of  2-D  for  the  2-pendulum,  there  are  two  distinct 
characteristic  angles  that  have  to  be  defined.  In  Figure  109,  the  angle  defined  by 
the  two  green  lines  will  be  called  the  forward  principal  characteristic  angle,  cpi, 
while  the  angle  defined  by  the  two  turquoise  lines  is  the  side  principal 
characteristic  angle,  cp2. 

1.  String  Tension  Investigation  of  the  Driven  4-Pendulum 

The  4-pendulum  is  assumed  to  be  subjected  to  a  forward,  lateral,  and 

vertical  support  motion  characterized  by  X ,  Z ,  and  Y  accelerations, 
respectively.  To  formulate  the  tension  equations  for  each  individual  string,  the 
superposition  rule  will  be  used.  Let  us  name  the  strings:  string  1,  the  one  closest 
to  the  global  system  of  axis;  and  strings  2-4,  counting  clockwise  looking  the 
pendulum  from  the  top  (Figure  109). 
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Figure  1 10.  4-pendulum  front  view  tension  breakdown 


Initially,  it  will  be  assumed  that  X  =  0 .  Thus,  the  analysis  will  be  similar  to 
the  one  completed  for  the  2-pendulum.  In  the  figure  above,  the  side  view  of  the 
4-pendulum  shows  that  the  strings,  represented  as  dashed  turquoise  lines,  lay  on 
the  same  plane  with  the  turquoise  non  vertical  lines  of  Figure  109.  If  one  of  these 
dashed  turquoise  lines  substituted  each  side  pair  of  strings,  the  principal  tension 
Tp  in  each  of  those  will  be  in  accordance  with  equations(2.62). 
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(3.1) 


Tpn 

Tp22 


mAY  +  g)  {  Z  . 
2  cos^2  sin^2 


rn^Y  +  g) 
2  cos^2 
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sin  ^ 
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Figure  111.  4-pendulum  side  view  tension  breakdown 


Looking  the  same  model  from  the  left-side  view,  Tp2\  can  be  analyzed 
into  two  components,  Tp2\\  and  Tp2\2,  since  they  correspond  to  string  1  and  2, 

respectively.  Let  us  name  the  angle  that  each  of  those  tension  components  forms 
with  the  resultant  tension,  (ppartiai2 ■  These  two  components  are  equal  and  their 
magnitude  is 


7>2i  1  =  7)321 2  =  Tp^1 - 

2  COS  ^partial! 


(3.2) 
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Then,  combining  equations  (3.1 )  and  (3.2), 

m 


Tp2\l  =  Tp2i2 


4  cos  (f) partjai2  cosrA  sin  <f>2 

Performing  the  same  analysis  for  strings  3  and  4  we  yield  to 


(3.3) 


Tp223=Tp224  ■ 


m 


( 


a  +  g) 


-) 


(3.4) 


4  cos  $ pariiai2  cos^2  sin^2 

If  we  rename  the  tensions  using  simpler  symbology  to  T1,  T2,  T3,  and  T4, 

respectively,  and  perform  the  same  methodology  for  the  case  that  Z  =0,  the 
yielding  results  are 

m 


T\  =  T4  = 


4cos^ 


partial  1 


cos^  sin^ 


T2  =  T3  = 


m 


4cos0 


( 


(X  +  g)  X 


(3.5) 


partial 1 


cos^  sin^ 


) 


Introducing  three  new  and  more  generalized  variables  H,  as  the  vertical 
position  of  the  mass  measured  along  the  negative  Y-axis,  Lfront,  as  the  frontal 
distance  between  the  support  points,  and  Lside,  as  the  side  view  distance 
between  the  support  points,  the  principal  characteristic  angles  can  be  expressed 
as  follows: 


,  ,  .Lside, 

<k  =  arctan(— ^— ) 

(/)2=avctML^lt) 


(3.6) 


Additionally,  let  us  define  Lpi  and  Lp2  ,  the  lengths  of  the  pseudo-strings 


of  the  principal  tensions  (green  and  turquoise  non-vertical  lines  in  Figure  101), 
and  / ,  the  length  of  each  string.  So, 
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(3.7) 


Therefore,  combining  (3.3),  (3.4),  (3.5)  and  (3.8),  the  tensions  exerted  on 
the  four  strings  for  every  combination  of  forward,  lateral  and  vertical 
accelerations  are 

n  = 


m,(Y  +  g)i  2X  2 Z  J^2  ,  fLJront^  ,  ( Lside  ^ 

A  v  ,,  "l_  J- _  •  |£ _ 77  "*"v  o  /  "*"v  O  / 


Lside  Lfront 


_  m  ((7  +  g)  2X  [  2Z  )  1^-2  !  (Lfront ~2  ,  ^ Lside ~2 
4  H  Lside  Lfront  V  2  2 
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This  set  of  equations  can  be  used,  as  already  done  for  the  corresponding 
one  of  the  2-pendulum,  to  predict  the  stability  of  the  4-pendulum,  i.e.,  at  which 
point  that  at  least  one  of  the  string  tensions  will  be  equal  to  zero. 

2.  Tension  Analysis  Validation 

Aiming  to  validate  the  above  set  of  equations,  the  VN  4-D  model  was 
used.  The  setup  of  the  model  was  complicated  due  to  the  required  level  of 
dimensional  precision.  Finally,  the  method  that  provided  the  most  reasonable 
results  (absence  of  nonlinearities  in  tension)  was  to  combine  two  separate  runs 
of  the  software  —  one  with  no  forward  acceleration  and  one  with  no  lateral 
acceleration  -  and  then  superimpose  the  results.  In  this  specific  model  /77=1kgr, 
the  frontal  distance  between  the  support  points  is  Lfront=  2m;  the  side  distance  is 
Ls/'c/e=1m;  and  the  vertical  position  of  the  mass  is  H=  2  m.  Finally,  the  applied 

accelerations  are  X  =  2?(m/sec2),  Z  =  4?(w/sec2),  and  7  =  0  (t  in  seconds). 


Figure  112.  Strings  1  and  2  tension  validation 
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Figure  113.  Strings  3  and  4  tension  validation 


B.  4-PENDULUM  DESIGN  OPTIMIZATION 

The  procedure  of  design  optimization  for  the  4-pendulum  requires  that  one 
defines  the  objective  function,  the  design  variables,  and  the  constraints  --  just  like 
it  was  done  for  the  2-pendulum.  Once  again,  the  objective  function  is  to  minimize 
the  tension  that  will  be  exerted  on  a  string  for  a  given  set  of  accelerations.  The 
design  variables  are  the  distances  H,  Lfront  and  Lside. 

According  to  (3.9)  set  of  equations,  the  Lfront ,  Lside ,  and  H  substitute 
<j\  and  (f)2  as  the  design  variables.  Nevertheless,  we  can  constrain  the  value  of  H 

with  an  equality  constraint,  instead  of  an  inequality  one,  to  be  able  to  obtain  a  2- 
D  visualization  of  the  objective  function. 

In  the  following  example,  the  known  values  are  H=  1m, 

Lfr0nt min  =  Lsidemm  =  ’  and 
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Figure  1 14.  Maximum  exerted  tension  for  arbitrary  accelerations 
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Figure  115.  Maximum  tensions  exerted  for  arbitrary  accelerations 


Figure  1 1 6.  Visualization  of  designs  that  fulfill  the  stability  constraint 
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Figure  117.  Final  visualization  of  maximum  tensions  for  arbitrary  accelerations 


In  Figure  117,  we  can  visualize  the  acceptable  Lfront  and  Lside 
combinations  and  at  the  same  time  find  the  optimum  solution  to  our  problem. 
Clearly  the  optimum  solution  is  Lfront=  Lside= 2m.  This  solution  can  also  be 
verified  using  the  fmincon  optimization  Matlab  function,  which  also  yields  the 
same  result. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  SUMMARY  AND  CONCLUSIONS 

The  following  table  summarizes  the  basic  equations  developed  throughout 


this  thesis. 


A/A 

Description 

Equation 

Equation 

No 

1 

Free 

oscillation 

of 

2- 

pendulum 

6  +  y  sin  9  cos  a  +  y  cos  6  sin  a  sign  9  =  0 

subject  to 

0+i  =  6_t  ,y(sin(90  -  2a))2  +  (r  cos(90  -  2a))2 

(2.18) 

and 

(2.55) 

2 

SCPA 

method 

prediction 

f  2  P  ) 

#  =  Acos  con  +  — - —  t 

V  71 M)  J 

2  g 

COq  =  yC0S<2 

P  =  ysincr 

(2.32) 

and 

(2.24) 

3 

Energy 

method 

prediction 

d  +  F  ^m6  =  0 

F  =  $  mg  sin (a  +  0)d6 1  mg  sin  Od 6 

(2.39) 

and 

(2.38) 

4 

Driven 

2- 

pendulum 

a  Ag  +  Y)  X  ■  ,  ■  n 

u  +  {^— — y_cosa  +  __  smcr)smfc'  + 

(( g  +  Y )  sina_^_  cosaysign0  cos  q  -  o 

(2.60) 
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Stability 
criterion  for 

driven 

2- 

pendulum 


X 


tan  a 


-Y>g 


(2.63) 


String 
tensions 
for  driven 

2- 

pendulum 
while  in 
equilibrium 
position 


n_mAY  +  g)  [  X  . 

2  v  cos  a  sin  a 
T2_m(jY  +  g)  X  ^ 
2  cos  a  sin  a 


(2.62) 


String 
tensions 
for  driven 

4- 

pendulum 
while  in 
equilibrium 
position 


2  Z 


n_ffl((7  +  g)|  2X 

4  H  Lside  Lfront 


2  Z 


)  \H2  +(Lfr°nt^2  ,  , Lside 


-)•  +( - ) 

2  y  2 


T2-m((Y+g)  2X 

4  H  Lside  Lfront 

T2~mAX+g)  2X  2  Z 

4  H  Lside  Lfront 

W-mAY+g)  |  2X  2Z 

4  H  Lside  Lfront 


)  I H2  +  {Lfwntf  +  ^ Lside y 


)  ifj2  +  (Mront\2.  ,  tLside ^2 


-)  +(-^r- ) 


)  l/f 2  +  (Lfront ^2  ,  (Lside ^ 


-)  +(-^— ) 


(3.9) 


Table  2.  Equations  summary 


Based  on  the  analysis  performed  throughout  this  thesis,  including  both  the 
analytical  formulation  of  the  mathematical  equations  and  their  validation 
performed  with  the  aid  of  VN  4-D,  we  can  draw  the  following  conclusions: 
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•  The  2-pendulum’s  configuration  of  two  strings  per  mass 
introduces  two  distinct  nonlinear  oscillation  characteristics.  The  first  one  results 
from  the  discontinuity  in  its  stiffness  around  the  equilibrium  position,  i.e.  an 
instantaneous  jump  in  the  value  of  the  restoring  moment  around  zero.  Four 
different  methods,  were  used  to  incorporate  this  discontinuity  in  analytical 
equation  formulation:  Numerical  solution  with  ode4  5,  exact  solution  for  a  time 
increment,  SCPA  method  and  Energy  method. 

•  The  second  nonlinear  oscillation  characteristic  is  related  to 
the  energy  dissipation  that  takes  place  at  the  switch  point.  The  dynamics  that 
govern  this  behavior  are  correlated  to  the  dynamics  of  a  standard  collision  of  a 
mass  on  a  rigid  wall,  with  specific  velocity  direction  constraints  before  and  after 
the  collision.  The  idealized  model  is  based  on  an  instantaneous  impact  and  total 
energy  loss  related  to  the  mass’  velocity  component  parallel  to  the  taught  string’s 
direction.  This  assumption  is  validated  by  the  VN  4-D  simulation. 

•  Another  characteristic  that  was  analyzed  was  the  stability  of 
the  2-string  pendulum  in  forced  oscillation.  This  string  configuration  results  in  the 
absence  of  oscillatory  behavior  if  the  mass  is  initially  at  the  equilibrium  position 
and  the  driving  horizontal  and  vertical  accelerations  do  not  exceed  values 
determined  exclusively  by  the  characteristic  angle  a.  This  observation  allows  the 
rational  design  of  support  systems  under  given  support  motion  characteristics. 

•  The  general  behavior  of  the  driven  2-pendulum  was 
investigated  including  the  case  of  induced  jumps.  Specific  examples 
demonstrated  the  highly  non-nonlinear  behavior  of  the  2-pendulum  when 
extreme  acceleration  values  of  the  support  cause  the  strings  to  go  slack. 

•  Finally,  the  analysis  of  the  2-pendulum  was  extended  on  the 
4-pendulum  configuration,  which  is  a  more  common  support  design  in  practice.  A 
typical  optimization  study  was  presented  in  order  to  visualize  the  practical  design 
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advantages  that  such  an  analysis  offers.  It  was  shown  that  it  is  possible  to 
reduce  both  motions  and  cable  tensions  by  appropriate  selections  of  string 
angles. 

B.  RECOMMENDATIONS 

Further  development  of  the  analysis  presented  in  this  work  could  be  in  the 
application  of  angular  support  motion,  both  for  the  case  of  forced  oscillations  of 
the  2-pendulum  and  the  trade  off  analysis.  This  work  would  provide  a  more 
realistic  representation  of  the  support  motion  analogous  to  the  motion  of  the  ship. 
Further  extensions  to  this  work  by  considering  the  effects  of  string  elasticity  are 
also  recommended. 
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APPENDIX: 


COMPUTER  PROGRAMS 


%Program  1 .m 

O, _ 

o 

o, _ 

o 

%Computing  the  response  of  free  2-pendulum  using  three  distinct 
methods:  Numerical  solution,  SCPA  method  and  Energy  method 

O, _ 

o 

o, _ 

o 


clear  all 
clc 

close  all 

global  a  omega2  A  omega  P  maxyl  Thetao  t  increased  restoring 
THETAo=pi/10  %initial  angular  displacement 

THETAo  dot=0  %initial  angular  velocity 


a=input ( ' Give  the  characteristic  angle  a  ' ) ; %2-pendulum' s 
characteristic 


g=9. 81; 
11=1; 


1=11/ cos (a) 
omega2=g/l; 

omega2o=omega2*cos (a)  ; 
tf=3  ; 


%angle 

%acceleration  of  gravity 

%equivalent  string  length  which  equals  to  the 
%simple  pendulum  string  length  or  the 
%perpendicular  distance  of  the  mass  from 
%the  plane  of  suspension 
%string  length  of  complex  pendulum 
%square  of  the  natural  frequency  of  the 
%  standard  pendulum 

%square  of  the  pseudo  natural  frequency  of  the 

%2-pendulum 

%final  time 


O, 

o 


%Numerical  solution  of  the  differential  equation  of  motion 


options  =  odeset ( ' RelTol ' , le-6, ' AbsTol ' , [ le-6  le-6]); 

[t, x] =ode45 ( ' Ffull ' , [0:0.01: tf ] , [THETAo, THETAo_dot] , options) ; 
[  t ,  x  (  :  ,  1 )  ]  ; 

XI  =x; 


figure  (1) 


hold  on 

title ( ' \alpha=\pi/6 ' ) 

plot (t, x ( : , 1) , 'g') 
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maxyl=max (x ( : , 1 ) ) ; 
maxy2=max (x ( : , 2) ) ; 


O, 

o 


%SCPA  method 


A=THETAo; 

B=sqrt (2*g* (11-cos (THETAo+a) ) / 1 )  ; 

t=0 : 0 . 01 : tf ; 

[xap, yap, Tapprox, Uo, Texact, tl , Ain] =SCPA (t, a, omega2 , A, omega, P, maxyl ) ; 


plot(t,xap, '--r') 


Xlap= [ ] ; 

Ylap= [ ] ; 

for  i=l : Tapprox* 1 00 ; 
t= ( i — 1 ) *0.01; 

[xap,  yap,  Tapprox,  Uo, Texact, tl , Ain] =SCPA (t, a, omega2 , A, omega, P, maxyl ) ; 

Xlap ( i ) =xap; 

Ylap ( i ) =yap; 
end 


%Energy  method 

o, _ 

o 


o, 

o 


syms  x  Theta  phi  increasing  factor  square  nat  freq  phi  t  Initial  theta 

Initial  theta  dot 

y=int ( ' sin (x) ',x,0, Theta) ; 

z=int ( ' sin (phi+x) ' , x, 0, Theta) ; 

z=subs ( z , phi , a) ; 

z=subs ( z , Theta, maxyl ) ; 

y=subs (y, Theta, maxyl ) ; 

increased_restoring=z/y; 


[t, y] =ode45 ('E',  [0:0.01: tf ]  ,  [THETAo, THETAo_dot] , options) ; 
[t, y  (  : ,  1)  ]  ; 
plot (t, y ( : , 1) , ' - . b ' ) 

Yl=y; 

Q, _ 

O 

o, _ 

o 
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t=0 : 0 . 01 : tf ; 


method ' ) 


hold  off 


xlabel ( ' t  (sec) ' ) 
ylabel ( ' \theta  (rad) ' ) 
grid  on 
box  on 

legend (' ODE45  ' , ' SCPA  method Energy 
axis([0  tf  -1.25*maxyl  0.6]) 


Tapproxl=Tapprox;  %  Approximate  Period  of  the  Slowly  Changing  method 
Texactl=Texact;  %Period  from  the  analytical  solution  calculating 

%  for  the  time  period  when  the  sign  of  theta  is 
%constant 


figure  (2) 

title ( ' \alpha=\pi/6 ' ) 
hold  on 

plot (XI ( 1 : Texactl* 1 00 , 1 ) . /THETAo, XI ( 1 : Texactl* 1 00 , 2 ) . / sqrt (omega2o) , ' g ' 
, Xlap . /THETAo, Ylap . / sqrt (omega2o) , ' --r ' , 

Y1 ( 1 : Tapproxl* 1 00 ,  1 )  . /THETAo, Y1 ( 1 : Tapproxl *  100 , 2 )  . / sqrt (omega2o) ,  ' - . b ' ) 
grid  on 
box  on 

xlabel ( ' \ theta / \thetao ' ) ; ylabel ( ' \ theta  1 ' / \ omega \o ' ) 
legend (' ODE45  ' , ' SCPA  method Energy  method') 
axis ([-1.5  1.5-1  1]  ) 

hold  off 
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%Subroutine  Ffull.m  of  Program  l.m  and  Program  2.m 

O, _ 

o 

o, _ 

o 

%Numerical  solution  of  the  differential  equation  of  motion  for 
%the  2-pendulum 

Q, _ 

o 

function  xp=Ffull (t, x) 
global  a  omega2 

xp=zeros (2,1)  ; 
xp ( 1 ) =x ( 2 ) ; 

xp(2)=-sin(x(l) ) *omega2*cos (a) -omega2*cos (x(l) ) * (sin(a) ) *sign(x(l) ) ; 


%Subroutine  E.m  of  Program  la.m  and  Program  lb.m 


O, 

o 


%Numerical  solution  of  the  differential  equation  of  motion  for 
%the  2-pendulum 

O, _ 

o 


function  yp=E(t,y) 

global  a  omega2  increased  restoring 

yp=zeros (2,1) ; 
yp (1) =y (2) ; 

yp  (2) =-sin (y (1) ) *omega2*increased_restoring; 


%Subroutine  SCPA.m  of  Program  l.m 


O, 

o 


%SCPA  method 


function 

[xap, yap, Tapprox, Uo, Texact, tl , Ain] =SCPA (t, a, omega2 , A, omega, P, maxyl ) ; 

global  a  omega2  A  omega  P  maxyl  Ain 

P=omega2*sin (a)  ; 

omega  =sqrt (omega2*cos (a) ) ; 

tl= (1/sqrt (omega2) ) *acos (1/ (1+ (A*omega2 *cos (a) ) / (omega2*sin (a) ) ) ) ; 
xap=A*cos ( (omega+2*P/ (pi*omega*A) ) *t) ; 

yap=A* (omega+2*P/ (pi*omega*A) ) *sin ( (omega+2*P/ (pi*omega*A) ) *t) ; 
Tapprox=2*pi/ (omega+2*P/ (pi*omega*A) ) ; 

Uo=-omega* (A+P/ (omega . A2 ) ) *sin (omega*tl ) ; 

Ain=-Uo/omega  - (P/ (omega. A2) ) *sin (omega*tl) ; 

Texact= (4/omega) *acos (1/ (l+A*omega*omega/P) ) ; 
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%Program  2 .m 


O, 

o 


%Computing  the  response  of  free  2-pendulum  using  three  distinct 
methods:  Numerical  solution,  SCPA  method  and  Energy  method 

O, _ 

o 


clear  all 
clc 

close  all 

global  a  omega2  A  omega  P  maxyl  Thetao  t  increased  restoring 
THETAo=0  %initial  angular  displacement 

THETAo  dot=1.22  %initial  angular  velocity 


a=input ( ' Give  the  characteristic  angle  a  ' ) ; %2-pendulum' s 
characteristic  angle 


g=9.81;  %acceleration  of  gravity 

11=1;  %equivalent  string  length  which  equals  to  the 

%simple  pendulum  string  length  or  the  %perpendicular  distance  of  the 
mass  from 

%the  plane  of  suspension 


1=11/ cos (a) 
omega2=g/ 1 ; 

omega2o=omega2 *cos (a) 

tf =3  ; 


%string  length  of  complex  pendulum 
%square  of  the  natural  frequency  of  the 
%  standard  pendulum 

%square  of  the  pseudo  natural  frequency  of  the 

%2-pendulum 

%final  time 


O, 

o 


%Numerical  solution  of  the  differential  equation  of  motion 


O, 

o 


options  =  odeset ( ' RelTol ' , le-6, ' AbsTol ' , [ le-6  le-6]); 

[t,  x] =ode45 ( ' Ffull '  ,  [0:0.01: tf ]  ,  [THETAo, THETAo_dot] , options) ; 
[t, x  (  :  ,  1) ]  ; 


%Predicting  the  maximum  displacement  based  on  the  initial  velocity 
%  using  energy  conservation  and  geometry 


Q. 

o 


so=l*cos (a) ; 
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h=THETAo_dot*THETAo_dot*l*l/ (2*g)  ; 
hl=ll-h; 

THETAo2=acos (hl/1) -a; 


figure  (1) 


hold  on 

title ( ' \alpha=\pi/6 ' ) 

plot  (t, x  ( : , 1) ,  ' g ' ) 


maxyl=max (x  ( : , 1 )  )  ; 


^SCPA  method 


A=THETAo2 ; 

t=0 : 0 . 01 : 3; 

[xap, Tapprox, Uo, Texact , tl , Ain] =SCPA  u ( t , a, omega2 , A, omega, P, maxyl ) ; 


plot(t,xap, ' — r ' ) 


%Energy  method 


syms  x  Theta  phi  increasing  factor  square  nat  freq  phi  t  Initial  theta 

Initial  theta  dot 

y=int ( ' sin (x) ' , x, 0 , Theta) ; 

z=int ( ' sin (phi+x) ' , x, 0 , Theta) ; 

z=subs ( z , phi ,  a)  ; 

z=subs ( z , Theta, maxyl )  ; 

y=subs (y, Theta, maxyl ) ; 

increased_restoring=z/y; 


[t, y] =ode45 ( ' E ' , [0:0.01: tf ] , [THETAo, THETAo_dot] , options) ; 
[t,  y  ( : ,  1)  ]  ; 
plot (t, y ( : , 1) , ' - . b ' ) 
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method ' ) 


hold  off 


xlabel ( ' t  (sec) ' ) 
ylabel ( ' \theta  (rad) ' ) 
grid  on 
box  on 

legend (' ODE45  ' , ' SCPA  method Energy 
axis([0  tf  -1.25*maxyl  3*maxyl]) 


%Subroutine  SCPA  U.m  of  Program  2.m 


O, 

o 


%SCPA  method  for  initial  velocity 


O, 

o 


function 

[xap, Tapprox, Uo, Texact, tl , Ain] =SCPA (t, a, omega2 , A, omega, P, maxyl ) ; 

global  a  omega2  A  omega  P  maxyl  Ain 

P=omega2*sin (a) ; 

omega  =sqrt (omega2*cos (a) ) ; 

tl= (1/ sqrt (omega2) ) *acos (1/ (1+ (A*omega2*cos (a) ) / (omega2*sin (a) ) ) ) ; 
xap=A*sin ( (omega+2*P/ (pi*omega*A) ) *t) ; 

Tapprox=2*pi/ (omega+2*P/ (pi*omega*A) ) ; 

Uo=-omega* (A+P/ (omega . A2 ) ) *sin (omega*tl ) ; 

Ain=-Uo/omega  - (P/ (omega. A2) ) *sin (omega*tl) ; 

Texact= ( 4 /omega) *acos ( 1 / ( 1+A* omega* omega/P) ) ; 
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%Program  3.m 


%Computing  the  response  of  the  driven  2-pendulum  for  a  given  and 
%  specific  set  of  support  motion 

O, _ 

o 


o, _ 

o 

9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-9-7\  rPrPT?'NTrr  TfiMS-  9-  S-9'S-9'S-9'S-9'S-S-S-9'S-9'S-9'S-9'S-9'S-9'S-&-S-S-S-9'S-9' 
oooooooooooooooooooooooooooooo  ±±J_iIM±J.\JiMoooooooooooooooooooooooooooooooo 

%this  m-file  is  valid  only  when  the  appropriate  and  specific  data_l.m 
%and  data  2.m  data  files  exist  in  the  same  folder  since  those  files 
%  contain  both  the  support  motion  info  and  the  corresponding  results 
%  from  Nastran  for  comparison  and  validation  purposes 


clc 

clear  all 
format  long 

fprintf ( ' make  sure  that  you  have  copied  the  data  provided  by  Nastran') 
datainput=input ( ' if  you  already  have  press  1  ') 

if  datainput==l 

fprintf (' good  to  go') 

else 

break 

end 


%call  the  data  file 

dat=input (' choose  the  data  file,  1  for  latteral  and  vertical  motion,  2 

for  latteral  motion') 

if  dat==l 

data_l 

else 

data_2 

end 


%response  before  the  stability  is  broken 

o, _ 

o 

o, _ 

o 


ti=0;  %initial  time  condition 

tf=length (data ( : , 2) ) /1000-0 . 001; %final  time  condition 
dt=0.001;  %time  step 
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nt=fix ( (tf-ti) /dt) 
a=pi/18 ; 
g=9. 81; 

11=1; 


1=11/ cos (a) ; 
omega2=g/ 1  ; 

omega2o=omega2*cos (a)  ; 

r=0  ; 
m=l  ; 


%number  of  time  steps 
%pendulation  characteristic  angle 
%acceleration  of  gravity 

%equivalent  string  length  which  equals  to  the 

%simple  pendulum  string  length  or  the 

%perpendicular  distance  of  the  mass  from  the 

%  plane  of  suspension 

%string  length  of  complex  pendulum 

%square  of  the  natural  frequency  of  the  simple 

%pendulum 

%square  of  the  pseudo  natural  frequency  of  the 
%2 -pendulum 

%restitution  coefficient 
%mass 


%initializing  the  acceleration,  velocity  and  displacement  vectors 

o, _ 

o 


acc= zeros ( 1 , nt+1 ) ; 
vel=zeros ( 1 , nt+1 ) ; 
disp= zeros ( 1 , nt+1 ) ; 
f or ce= zeros ( 1 , nt+1 ) ; 
Zabs=zeros (l,nt+l)  ; 
Yabs=zeros (l,nt+l)  ; 
Z=zeros ( 1 , nt+1 ) ; 
Y=zeros ( 1 , nt+1 ) ; 
Zs=zeros ( 1 , nt+1 ) ; 
Ys=zeros ( 1 , nt+1 ) ; 
acc_rest= zeros ( 1 ,  nt+1 ) 
%that  the  mass  doesn't 
Zsddot=zeros ( 1 ,  nt+1 )  ; 
%midpoint 

Ysddot=zeros ( 1 , nt+1 ) ; 
%midpoint 

Zsdot= zeros ( 1 , nt+1 ) ; 
Ysdot= zeros ( 1 , nt+1 ) ; 
Zsdot ( 1 , 1 ) =0 ; 

Ysdot ( 1 , 1 ) =0 ; 

Zsddot ( 1 , 1) =0 ; 

Ysddot (1 , 1) =0 ; 
Vzj=zeros ( 1 , nt+1 ) ; 

Vyj  = zeros ( 1 , nt+1 ) ; 


%mass  angular  acceleration 
%mass  angular  velocity 
%mass  angular  displacement 
%pseudo-restoring  force  on  mass 
%absolute  horizontal  position  of  mass 
%absolute  vertical  position  of  mass 
%global  horizontal  position  of  mass 
%global  vertical  position  of  mass 
%global  horizontal  position  of  support  midpoint 
%global  vertical  position  of  support  midpoint 
%pseudo-acceleration  corresponding  to  the  case 
break  stability 

%global  horizontal  acceleration  of  support 

%global  vertical  acceleration  of  support 

%global  horizontal  velocity  of  support  midpoint 
%global  vertical  velocity  of  support  midpoint 


%linear  horizontal  velocity  of  mass  during  jump 
%linear  vertical  velocity  of  mass  during  jump 


O, 

o 


%applying  the  data  input  in  the  support  midpoint  motion 


if  dat==l 

Zs=data ( : , 2) ' . /1000; 
Ys=data ( : , 3) ' . /1000; 
else 

Zs=data ( : , 2)  '  .  /1000; 
Ys=zeros ( 1 , nt+1 ) ; 
end 
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%initializing  the  tension,  support  points  position  and  string  length 
%vectors 


O, 

o 


Tright=zeros ( 1 , nt+1 )  ; 
Tie ft= zeros ( 1 ,  nt+1 )  ; 
Zright= zeros ( 1 , nt+1 ) ; 
Zlef t= zeros ( 1 , nt+1 ) ; 
Yright= zeros ( 1 , nt+1 ) ; 
Ylef t= zeros ( 1 , nt+1) ; 
Lright=zeros ( 1 , nt+1 ) ; 
Lief t= zeros ( 1 , nt+1 ) ; 


%right  string  tension 
%left  string  tension 

%horizontal  position  of  right  support  point 
%horizontal  position  of  left  support  point 
%vertical  position  of  right  support  point 
%vertical  position  of  left  support  point 
%distance  of  mass  from  right  support  point 
%distance  of  mass  from  left  support  point 


for  it=l:nt-l;  %calculating  the  support  displacements  and 

%accelerations  based  on  the  support  motion 
Zsdot (1,  it+1)  =  (Zs (1, it+1) -Zs  (1, it) )  /dt; 

Ysdot (1, it+1) = (Ys (1, it+1) -Ys (1, it) ) /dt; 

Zsddot (1, it+1)  =  (Zsdot (1, it+1) -Zsdot (1,  it) ) /dt; 

Ysddot (1, it+1) = (Ysdot (1, it+1) -Ysdot (1, it) ) /dt; 

end 

vel(l,l)=0;  %initial  angulat  velocity  condition 

disp(l,l)=0;  %initial  angular  displacement  condition 

Zabs(l,l)=-1* (sin (abs (disp (1, 1) ) +a) -sin(a) ) ; 

Yabs(l,l)=ll-1* (cos (abs (disp (1, 1)  )+a)  )  ; 

if  abs (max (Ysddot) ) ==0  &&  abs (max (Zsddot) ) ==0&&  abs (disp ( 1 , 1 ) ==0 ) 
%if  there  is  no  acceleration  applied  to  the  support  points 
%there  is  no  delay  for  breaking  the  stability 

T=1  ; 

else 


for  it=l:nt; 

acc_rest (1,  it)  =  (abs (Zsddot (1, it) ) ) /tan  (a)  -Ysddot (1, it) ; 

Tright ( 1 , it) = (m/2 ) * ( (Ysddot (1, it) +g) /cos (a)  + Zsddot ( 1 , it) /sin  (a) ) ; 
Tleft(l,it)  =  (m/2)*( (Ysddot ( 1 ,  it) +g) /cos (a)  -Zsddot (1,  it) /sin (a) )  ; 

Z right (l,it)=Zs(l,it)+l*tan(a)  ; 

Zleft  (1, it) =Zs (1, it) -l*tan (a)  ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

if  acc_rest ( 1 , it) <g  %if  the  combined  applied 

%acceleration  is  less  than  the  required  to  break  the  stability 

disp (1, it) =0;  %there  is  no  oscillation 

vel (1, it) =0; 
acc (1, it) =0; 

Zabs (1,  it) =-l* (sin (disp (1, 1) +a) -sin  (a) ) ; 

Y abs (l,it)=ll-l* (cos (abs (disp (1,1) )+a) ) ; 
Z(l,it)=Zabs(l,it)+Zs(l,it) ; 

Y  (1,  it)  =Yabs  (1,  it)  +Ys  (1,  it)  -11; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) .A2+(Y(l,it)- 
Yright (1, it) )  . A2)  ; 
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Lleft (1, it) =sqrt ( (Z(l,it)-Zleft(l,it) ) ,A2+(Y(l,it)~ 
Yleft (1, it) )  .  A2)  ; 

else 

T=it; 

Zabs  (1,  it)  =-l*  (sin  (disp  (1,  1)  +a)  -  sin  (a)  )  ; 

Yabs (1, it)  KL1-1* (cos (abs ( disp (1,1) )+a) ) ; 
Z(l,it)=Zabs(l,it)+Zs(l,it) ; 

Y (1,  it) =Yabs (1, it) +Ys (1, it) -11; 

Lright (1,  it) =sqrt ( (Z (1, it) -Zright  (1, it) )  ,A2+(Y(l,it)~ 
Yright (1, it) )  . A2)  ; 

Lleft (1,  it) =sqrt ( (Z (1, it) -Zleft (1, it) )  ,A2+(Y(l,it)~ 
Yleft (1, it) )  .  A2)  ; 

Vzl=Zsdot (1, T) ; 

Vyl=Ysdot (1, T) ; 

ACC=Zsddot (1, T) ; 

break 

end 


end 


end 


o, 

o 


%Step  after  the  stability  is  broken 


for  it=T:nt;  %loopl  -  initiating  calculations  after 

%stability  is  brocken 

if  (Tright ( 1 ,  it)  >0  |  |  Tleft (1, it) >0) %conditionl  -  oscillations  take 
%place  only  if  at  least  one  of  the  tensions  has  a  non-zero  value 
force ( 1 , it) =- 

(omega2+Ysddot (1, it) / 1 ) * (sin(a) ) *cos (disp (1, it) ) *sign (disp (1,  it) )  ; 
acc ( 1 , it) =force ( 1 , it) - 

(omega2+Ysddot  (1,  it)  /l)  *  sin  (disp  (1,  it)  )  *cos  (a)  +  (Zsddot  ( 1 ,  it)  )  *  (  (cos  (abs 
(disp  (1,  it)  )  +a)  )  )  /l; 
vel ( 1 ,  it+1 ) =vel (1, it) face (1, it) *dt; 

disp (1,  it+1) =disp (1, it) +vel (1, it+1) *dt; 

if  disp ( 1 , it+1 ) *disp ( 1 , it) >=0  %condition2  -  calculations 

%while  the  mass  oscillates  without  changing  strings 
disp (1,  it+1) =disp (1, it) +vel (1, it+1) *dt; 

Zabs(l,it+l)=-sign( (disp (1, it+1) ) ) *1* (sin(abs(disp(l,it+l) ) +a) - 
sin (a) ) ; 

Yabs (1, it+1) =11-1* (cos (abs (disp (1, it+1) )+a) ) ; 

Z (1,  it+1) =Zabs (1, it+1) +Zs (1, it+1) ; 

Y (1,  it+1) =Yabs (1, it+1) +Ys (1, it+1) -11; 

Zright (l,it)=Zs(l,it)+l*tan(a)  ; 

Zleft (1, it) =Zs (1, it) -1* tan (a) ; 

Yright (1, it) =Ys (1, it) ; 

Yleft(l,it)=Ys(l,it) ; 
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Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) .A2+(Y(l,it)- 
Yright (1, it) )  .  A2)  ; 

Lleft (1,  it) =sqrt ( (Z (1, it) -Zleft  (1, it) )  ,A2+(Y(l,it)- 
Yleft (1, it) )  .  A2)  ; 


else  %condition2  -  the  mass  reached 

%a  string  switch  point 
Tl=it; 

V=vel (1, Tl) ; 

Vaft=sign (V) *sqrt ( (V*sin (pi/2  -2*a) ) . A2+ (r*V*cos (pi/2  -2*a) ) . A2) ; 
vel ( 1 , Tl+1 ) =Vaft; 

disp (1, Tl  +  1) =disp (1,  it) +vel (1,  Tl  +  1) *dt; 

Zabs(l,it+l)=-sign( (disp ( 1 , it+1 ) ) ) *1* (sin(abs(disp(l,it+l) )  +a)  - 
sin (a) ) ; 

Yabs(l, it+1) =11-1* (cos (abs (disp ( 1 , it+1 ) )+a) ) ; 

Z (1,  it+1) =Zabs (1, it+1) +Zs (1, it+1) ; 

Y (1, it+1) =Yabs (1, it+1) +Ys (1, it+1) -11; 

Zright (l,it)=Zs(l,it)+l*tan(a)  ; 

Zleft (1, it) =Zs (1, it) -1* tan (a) ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) ,A2+(Y(l,it)- 
Yright (1, it) ) . A2) ; 

Lleft (1,  it) =sqrt ( (Z (1, it) -Zleft  (1, it) )  ,A2+(Y(l,it)- 
Yleft (1, it) )  .  A2)  ; 

end  %condition2 

if  disp ( 1 , it+1 ) >0  %condition3  -  calculating 

%the  strings  tension  depending  on  which  one  of  the  strings  is  taught 

Tright (1, it+1) =m* (vel (1, it+1) *vel (1, it+1) *l+Zsddot (1, it+1) *sin (abs (disp 
(1, it+1) ) +a) + (g+Ysddot (1, it+1) )* cos (abs (disp (1, it+1) )+a) ) ; 

Tleft (1, it+1) =0; 
if  Tright (1, it+1) >0 

Tright ( 1 ,  it+1 ) =T right ( 1 ,  it+1 )  ; 
else 

Tright (1 , it+1) =0  ; 

end 

else  %condition3 

Tleft (1, it+1) =m* (vel (1, it+1) *vel (1,  it+1) *1- 
Zsddot (1, it+1) *sin(abs (disp (1, it+1) )+a)+ (g+Ysddot (1, it+1) ) *cos (abs (disp 
(1, it+1) ) +a) ) ; 

Tright (1 , it+1) =0; 
if  Tlef t ( 1 , it+1 ) >0 

Tleft (1, it+1) =Tlef t ( 1 , it+1) ; 
else 

Tleft (1, it+1) =0; 

end 

end  %condition3 

else  %conditionl  -  the  first 

%jump  starts 
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T2=it 


V=vel ( 1 , T2 ) ; 

disp (1,  T2  +  1) =disp (1, T2) +vel (1, T2) *dt; 

Zabs(l,it)=-sign( (disp (1, it) ) ) *1* (sin (abs (disp (1, it) ) +a) -sin (a) ) ; 
Yabs(l,it+l)=ll-l* (cos (abs (disp (1,  it)  )+a)  )  ; 

Z  ( 1 , it) =Zabs (1 , it) +Zs (1 , it) ; 

Y  (1,  it)  =Yabs  (1,  it)  +Ys  (1,  it)  -11; 

Z right  (1,  it)  =Zs  (1,  it)  +l*tan  (a)  ; 

Zleft(l,it)=Zs(l,it)-l*tan(a)  ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) ,A2+(Y(l,it)- 
Yright (1, it) )  .  A2)  ; 

Lleft (1,  it) =sqrt ( (Z(l,it)-Zleft(l,it) )  ,A2+(Y(l,it)- 
Yleft (1, it) )  .  A2)  ; 

if  disp(l,T2)>0  %condition4  -  calculating 

%the  linear  velocity  of  the  mass  at  the  moment  where  the  jump  initiates 
Vmagn=abs ( l*vel ( 1 , T2 ) ) ; 

Vz2=-sign (V) *Vmagn*cos ( (disp ( 1 , T2 ) +a) ) +Zsdot ( 1 , T2-1 ) ; 

Vy2=sign (V) *Vmagn*sin ( (disp ( 1 , T2 ) +a) ) +Ysdot ( 1 , T2-1 ) ; 
break 

else  %condition4 

Vmagn=l*abs (vel ( 1 , T2 ) ) ; 

Vz2=-sign (V) *Vmagn*cos ( (abs (disp ( 1 , T2 ) )+a) )+Zsdot(l,T2-l) ; 
Vy2=-sign (V) *Vmagn*sin ( (abs (disp ( 1 , T2 ) ) +a) ) +Ysdot ( 1 ,  T2-1 ) ; 
break 

end  %condition4 

end  %conditionl 


end  %loopl 

%correcting  the  computation  rounding  errors 
LLright=Lright (1, T2) -0 . 004; 

LLlef t=Llef t ( 1 , T2 ) -0 . 004 ; 

Lright ( 1 , T2 ) =Lright ( 1 , T2 ) -0 . 004 ; 

Lleft (1, T2) =Lleft (1, T2) -0.004; 


while  it<nt 
j  ump 

af ter j  ump 
end 


time  =data  ( : , 1 )  ; 
if  dat==l 
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figure  (4) 

plot (time, disp, time, data ( :  ,  4 )  ,  ' — r  '  ) 
xlabel ( ' t  (sec) ' ) ; ylabel ( ' \ theta ( rad) ' ) ; 
axis([0  tf  -pi/2  pi/2]) 
grid  on 
figure  (5) 

subplot (2, 1, 1) ;plot (time, Zabs , time, data ( : , 1 0 ) ./1000, ' --r '  ) 

xlabel ('t  (sec) '); ylabel (' absolute  horizontal  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  1] ) 

subplot (2 , 1 , 2 ); plot (time, Yabs , time, data (:,  9)  . /1000  +  11  ,'--r'  ) 
xlabel ('t  (sec) '); ylabel (' absolute  vertical  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  1] ) 
figure  (6) 

subplot (2 , 1 , 1 ); plot (time, Z , time, data (:, 8 ). /1000  ,'--r'  ) 
xlabel ('t  (sec)  '); ylabel (' total  horizontal  position  (m)  ' ) ; 
grid  on 

axis (  [0  tf  -1  3] ) 

subplot (2 , 1 , 2 ); plot (time, Y, time, data (:,  7 ). /1000  ,'--r'  ) 
xlabel ('t  (sec) '); ylabel (' total  vertical  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  1] ) 
else 

figure  (4) 

plot (time, disp, time, data (:,3), ' --r ' ) 
xlabel ( ' t  (sec)  ' ) ; ylabel ( ' \theta (rad)  ' )  ; 
axis([0  tf  -pi/2  pi/2]) 
grid  on 
figure  (5) 

subplot (2, 1, 1) ;plot (time, Zabs , time, data (:,4) ,/1000-Zs', ' --r '  ) 

xlabel ('t  (sec) '); ylabel (' absolute  horizontal  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  1] ) 

subplot (2, 1,2);  plot (time, Yabs , time, data (:,5)./1000  +  ll-Ys'  ,'--r'  ) 
xlabel ('t  (sec) '); ylabel (' absolute  vertical  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  1] ) 
figure  (6) 

subplot (2 , 1 , 1 ); plot (time, Z , time, data (:, 4 ). /1000  ,'--r'  ) 
xlabel ('t  (sec) '); ylabel (' total  horizontal  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  5] ) 

subplot (2 , 1 , 2 ); plot (time, Y, time, data (:,  5 ). /1000  ,'--r'  ) 
xlabel ('t  (sec) '); ylabel (' total  vertical  position  (m) ' ) ; 
grid  on 

axis (  [0  tf  -1  0] ) 
end 

figure  (7) 

subplot (2,1,1) ;plot (time, Lright  ) 

xlabel ('t  (sec) '); ylabel (' absolute  horizontal  position  (m) ' ) ; 
grid  on 


126 


axis ( [ 0  tf  0  2 ] ) 

subplot (2,1,2) ;plot (time,  Lleft  ) 

xlabel('t  (sec) '); ylabel (' absolute  vertical  position  (m) ' ) ; 
grid  on 

axis (  [ 0  tf  0  2 ] ) 


figure  (8) 

subplot (2, 1, 1) ;plot (time, Zs  ) 

xlabel('t  (sec)  '); ylabel ( '  horizontal  support  position  (m)  ' ) ; 
grid  on 

axis (  [ 0  tf  0  5 ] ) 

subplot (2, 1, 2) ;plot (time, Ys  ) 

xlabel('t  ( sec)  '); ylabel (' vertical  support  position  (m)  ' ) ; 
grid  on 

axis (  [0  tf  -1  2] ) 


figure  (10) 

plot (time, Tright,  time,  Tleft) 

xlabel ( ' t  ( sec)  ' )  ; 

ylabel (' tension  (N) ') 

legend (' right  string left  string') 

axis  tight 


%jump.m  subroutine  -  part  of  the  Program  3.m  code 

O, _ 

o 


%Computing  the  response  of  the  driven  2-pendulum  when  a  jump  is 
initiated 


O, 

o 


o, 

o 


%collecting  data  from  the  main  file 


Lright ( 1 , T2 ) =LLright; 
Lleft (1, T2) =LLleft; 
Vyj ( 1 , T2 ) =Vy2 ; 

Vzj  (1,  T2) =Vz2; 
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loopl  -  initiating 


for  it=T2+l:nt; 

%calculations 

if  (Lleft ( 1 , it-1 ) >1  I  I  Lright (1, it-1) >1)  %conditionl  -  if  the 

%length  of  at  least  one  of  the  strings  is  equal  to  1,  the  the  jump  is 
%terminated 

T3=it; 

Z right (1,  it) =Zs (1, it) +l*tan  (a) ; 

Zleft(l,it)=Zs(l,it)-l*tan(a) ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

V  z  j  ( 1 ,  i  t )  =  V  z  2  ; 

Vyj (1, it) =Vyj (1 , it-1 ) - (g) *dt; 

Z (1,  it) =Z (1, it-1) +Vz2*dt; 

Y (1, it) =Y (1, it-1) +Vyj (1, it) *dt; 

Zabs(l,it)=Z(l,it)-Zs(l,it) ; 

Yabs  (1,  it)  =11 +Y  (1,  it)  -  Ys  (1,  it)  ; 

vel (1, it) =0; 

disp (1, it) =0; 

V  z  j  1  = V  z  j  ( 1 ,  i  t )  ; 

Vyj l=Vyj (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) ,A2+(Y(l,it)~ 

Yright (1, it) ) .A2) -0.004; 

Lleft (1, it) =sqrt ( (Z(l,it)-Zleft(l,it) ) ,A2+(Y(l,it)~ 

Yleft  (1, it) )  . A2) -0.004; 

Vmagn j  =sqrt (Vz j l*Vz j 1+Vyj 1+Vyj 1 ) ; 

break 

else 

Vz j ( 1 , it) =Vz2 ;  %conditionl  -  the  mass 

%follows  a  projectile  path 

Vyj  ( 1 , it) =Vyj  ( 1 ,  it-1 ) - (g) *dt;  %during  the  jump 

Z (1, it) =Z ( 1 , it-1) + V z  2  *  dt ; 

Y (1, it) =Y (1, it-1) +Vyj ( 1 , it) *dt; 

Zabs(l,it)=Z(l,it)-Zs(l,it) ; 

Yabs  (1,  it)  =11 +Y  (1,  it)  -  Ys  (1,  it)  ; 

Zright (l,it)=Zs(l,it)+l*tan(a) ; 

Zleft(l,it)=Zs(l,it)-l*tan(a)  ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) .A2+(Y(l,it)- 
Yright (1, it) )  .A2) -0.004  ; 

Lleft (1, it) =sqrt ( (Z (1, it) -Zleft (1, it) ) .A2+(Y(l,it)- 
Yleft  (1, it) )  .A2) -0.004 ; 

vel (1, it) =0; 

disp (1, it) =0; 

end  %conditionl 

end  %loopl 


if  Lright ( 1 , T3 ) >Lleft ( 1 ,  T3 ) 

%condition2  -  calculating  the  angular  displacement 
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phi=atan ( ( Z right (1,T3) -Z (1,T3) ) / ( Yright ( 1 , T3 ) -Y (1,T3) ) )  ;  %and 
%the  velocity  of  the  mass  after  the 

alpha=atan (Vyj 1/Vz j 1 ) ;  %jump 

%is  terminated  depending  on  which  string  is  taught 

vel (1,T3)=(1/1) * (-Vzjl*cos (phi) +Vyj 1* sin (phi) ) ; 
disp ( 1 , T3 ) =phi-a; 

Tright (1, T3) =m* (vel (1,T3) *vel (1, T3) *l+Zsddot (l,T3)*sin(abs(disp(l,T3) )+ 
a) + (g+Ysddot (1, T3) ) *cos (abs (disp (1, T3) ) +a) ) ; 

Tleft (1,  T3) =0; 

if  Tright ( 1 , it+1 ) >0  %condition3  -  ensuring 

%that  the  string  doesn't  support  compressive  loads 
Tright ( 1 , it+1 ) =Tright ( 1 , it+1 ) ; 

else 

Tright ( 1 , it+1 ) =0 ; 

end 

else  %condition2 

phi=atan ( (Z (1, T3) -Zleft (1,T3) ) / (Yleft (1,T3) -Y (1, T3) ) ) ; 
phit=-phi ; 

alpha=atan (Vz j 1/Vyj 1 ) ; 
alpha*180/pi; 

vel  (1, T3)  =  (1/1) * (-Vzjl*cos (phi) -Vyj 1* sin (phi)  )  ; 
disp ( 1 , T3 ) =- (phi -a) ; 

Tleft (1, T3) =m* (vel (1, T3) *vel (1,T3) *l+Zsddot (1, T3) *sin(abs (disp (1, T3) )+a 
)  +  (g+Ysddot (1, T3) ) *cos (abs (disp  (1, T3) ) +a)  )  ; 

Tright (1, T3) =0; 

if  Tlef t ( 1 , it+1 ) >0  %condition3 

Tleft (1, it+1) =Tlef t ( 1 , it+1) ; 

else 

Tleft (1, it+1) =0; 

end 

end  %condition2 
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%afterjump.m  subroutine  -  part  of  the  Program  3.m  code 


%Computing  the  response  of  the  driven  2-pendulum  after  a  jump  is 
%terminated 


O, 

o 


for  it=T3:nt;  %loopl  -  initiating  calculations  after  stability  is 
%broken 

if  (Tright ( 1 ,  it) >0  |  |  Tleft (1,  it)  >0)  %conditionl  -  oscillations  take 
%place  only  if  at  least  one  of  the  tensions  has  a  non-zero  value 

force ( 1 , it) =- 

(omega2+Ysddot (1, it) /l) * (sin(a) ) *cos (disp (1, it) ) *sign (disp (1, it) ) ; 
acc ( 1 , it) =force (1, it)  - 

(omega2+Ysddot (1, it) /l) *sin (disp (1,  it) ) *cos (a)  +  (Zsddot (l,it))*((cos(abs 
(disp  (1,  it)  )  +a)  )  )  /l; 
vel ( 1 ,  it+1 ) =vel (1, it) +acc (1, it) *dt; 

disp (1,  it+1) =disp (1, it) +vel (1, it+1) *dt; 

if  disp ( 1 , it+1 ) *disp ( 1 , it) >=0  %condition2  -  calculations 

%while  the  mass  oscillates  without  changing  strings 

disp (1,  it+1) =disp (1, it) +vel (1, it+1) *dt; 

Zabs(l,it+l)=-sign( (disp (1, it+1) ) ) *1* (sin(abs(disp(l,it+l) ) +a) - 
sin (a) ) ; 

Yabs (1, it+1) =11-1* (cos (abs (disp (1, it+1) )+a) ) ; 

Z (1,  it+1) =Zabs (1, it+1) +Zs (1, it+1) ; 

Y (1, it+1) =Yabs (1, it+1) +Ys (1, it+1) -11; 

Z right  (1,  it)  =Zs  (1,  it)  +l*tan  (a)  ; 

Zleft(l,it)=Zs(l,it)-l*tan(a)  ; 

Yright (l,it)=Ys(l,it); 

Ylef t (1, it) =Ys (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) ,A2+(Y(l,it)- 
Yright (1, it) ) .  A2)  ; 

Lleft (1,  it) =sqrt ( (Z(l,it)-Zleft(l,it) )  ,A2+(Y(l,it)~ 

Yleft (1, it) )  .  A2)  ; 


else  %condition2  -  the  mass  reached 

%a  string  switch  point 

T2=it; 

V=vel (1, T2) ; 

Vaft=sign (V) *sqrt ( (V*sin (pi/2  -2*a) ) . A2+ (r*V*cos (pi/2  -2*a) ) . A2) ; 
vel ( 1 , T2+1 ) =Vaft; 

disp ( 1 ,  T2  +  1 ) =disp ( 1 , T2 ) +vel ( 1 , T2  +  1 ) *dt; 

Zabs(l,it+l)=-sign( (disp (1, it+1) ) )*1* (sin(abs(disp(l,it+l) ) +a) - 
sin (a) )  ; 
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Yabs(l, it+1) =11-1* (cos (abs (disp ( 1 , it+1 ) )+a) ) ; 

Z (1,  it+1)  =Zabs (1, it+1) +Zs (1, it+1) ; 

Y (1, it+1) =Yabs (1, it+1) +Ys (1, it+1) -11; 

Z right (l,it)=Zs(l,it)+l*tan(a) ; 

Zleft  (1, it) =Zs (1, it) -l*tan (a)  ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) ,A2+(Y(l,it)- 
Yright (1, it) ) .  A2)  ; 

Lleft (1,  it) =sqrt ( (Z (1, it) -Zleft  (1, it) )  ,A2+(Y(l,it)- 
Yleft (1, it) )  .  A2)  ; 


end  %condition2 

if  disp ( 1 , it+1 ) >0  %condition3  -  calculating 

the  strings  tension  depending  on  which  one  of  the  strings  is  taught 

Tright ( 1 , it+1 ) =m* (vel ( 1 , it+1 ) *vel ( 1 , it+1 ) *l+Zsddot (1, it+1) *sin (abs (disp 
(1, it+1) ) +a)  +  (g+Ysddot (1, it+1) )*cos (abs (disp (1, it+1) )+a) )  ; 

Tleft (1, it+1) =  0 ; 

if  Tright ( 1 , it+1 ) >0  %condition4  -  ensuring  that 

%the  string  doesn't  support  compressive  loads 

Tright ( 1 , it+1 ) =Tright ( 1 , it+1 ) ; 
else 

Tright (1 , it+1) =0 ; 

end 

else  %condition3 

Tleft (1,  it+1) =m* (vel (1, it+1) *vel (1, it+1) *1- 
Zsddot (l,it+l)*sin(abs(disp(l,it+l))+a)+ (g+Ysddot ( 1 , it+1 ) ) *cos (abs (disp 
(1, it+1) ) +a) ) ; 

Tright (1 , it+1) =0 ; 

if  Tlef t ( 1 , it+1 ) >0  %condition4  -  ensuring  that 

%the  string  doesn't  support  compressive  loads 

Tleft (1, it+1) =Tlef t ( 1 , it+1) ; 

else 

Tleft (1, it+1) =0 ; 

end 

end  %condition3 

else  %conditionl 

T2=it 

V=vel ( 1 , T2 ) ; 

disp  ( 1 ,  T2  +  1 ) =disp ( 1 , it) +vel ( 1 , T2 ) *dt; 

Zabs(l,it)=-sign( (disp (1, it) ) ) *1* (sin (abs (disp (1, it) ) +a) -sin(a) ) ; 
Yabs(l, it+1) =11-1* (cos (abs (disp ( 1 , it) )+a) ) ; 
Z(l,it)=Zabs(l,it)+Zs(l,it) ; 
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Y (1,  it)  =Yabs (1, it) +Ys (1, it) -11; 

Z right  (1, it) =Zs (1, it) +l*tan (a)  ; 

Zleft  (1, it) =Zs (1, it) -l*tan (a) ; 

Yright (1, it) =Ys (1, it) ; 

Yleft (1, it) =Ys (1, it) ; 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (1, it) ) ,A2+(Y(l,it)- 
Yright (1, it) )  .  A2)  ; 

Lleft (1,  it) =sqrt ( (Z (1, it) -Zleft  (1, it) )  .A2+(Y(l,it)- 
Yleft (1, it) )  .  A2)  ; 

if  disp ( 1 , T2+1 ) >0  %condition5  -  calculating 

%the  velocity  when  ont  of  the  strings  becomes  taught 
Vmagn=abs ( l*vel ( 1 , T2 ) ) ; 

Vz2=-sign (V) *Vmagn*cos ( (disp ( 1 ,  T2 ) +a) ) +Zsdot ( 1 ,  T2-1 )  ; 

Vy2=sign (V) *Vmagn*sin ( (disp ( 1 , T2 ) +a) ) +Ysdot ( 1 , T2-1 ) ; 
break 

else  %condition5 

Vmagn=abs ( l*vel ( 1 , T2 ) ) ; 

Vz2=-sign (V) *Vmagn*cos ( (abs (disp (1, T2) ) +a) ) +Zsdot (1, T2-1) ; 
Vy2=-sign (V) *Vmagn*sin ( (abs (disp ( 1 , T2 ) ) +a) ) tYsdot ( 1 , T2-1 ) ; 
break 

end  %condition5 

end  %conditionl 


end  %loopl 

vel ( 1 ,  T2  )  ; 

Vz2  ; 

Vy2  ; 

LLright=Lright (1, T2) -0 . 004; 

LLlef t=Llef t ( 1 ,  T2 ) -0.004; 
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